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Preface 


This  report,  presents  a  code  that  calculates  the  neutron 
and  gamma  flueiices  (or  doses)  that  restit  from  a  nuclear 
weapon  explosion  between  7  and  100  kilometers  altitude  in  an 
exponential  atmosphere.  The  code  also  determines  if  air¬ 
craft  located  in  the  vicinity  survive  the  neutron  and  gamma 
effects.  The  code  is  listed  in  Appendix  B  and  the  instruc¬ 
tions  for  using  the  code  are  in  Chapter  III. 

The  purpose  of  developing  this  code  was  to  provide  a 
subroutine  for  gamma  and  neutron  effects  to  a  general  nuclear 
effects  survivability/vulnerability  code.  I  was  successful 
in  developing  such  a  code;  however,  the  time  required  to 
run  it  does  limit  its  applications.  The  time  required  is  about 
20  minutes  central  processor  time  and  an  hour  of  input/output 
time . 


Many  pev,nle,  in  the  course  of  my  efforts  to  complete 
this  code,  provided  assistance.  These  people  were  from  the 
Air  Force  Weapons  Laboratory  (AFWL) ,  Oak  Ridge  National 
Laboratory,  and  the  Air  Force  Institute  of  Technology  (AFIT) . 
To  all  of  you  I  offer  my  thanks.  I  would  like  to  single  out 
a  few  for  special  thanks,  because  without  their  help,  I  could 
not  have  completed  this  task.  First  is  Mr.  Harry  Murphy  of 
AFWL  for  his  talks  with  me  on  the  code,  SMAUG.  I  borrowed 
heavily  from  this  code  to  develop  mine.  Mr.  Robert  Roussin 
of  the  Radiation  Shielding  Information  Center  of  Oak  Ridge 
furnished  the  cross  section  data  I  used  in  the  code.  Next 
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are  two  members  of  the  AF1T  Mathematics  Department,  Major 
John  Jones  and  Dr.  Wilhelm  Ericksen.  Major  Jones  and  I  spent 
considerable  time  investigating  method:;  of  directly  solving 
sparse  matrix  equations.  Dr.  Ericksen  provided  the  basis 
for  the  nonortbogonal  coordinate  system  I  used.  Dr.  Donn 
Shankland  of  the  AFIT  Physics  Department  provided  f  irther 
assistance  on  this  coordinate  system  and  removed  the  last 
problem  area  standing  in  the  way  of  success.  I  can  not,  of 
course,  forget  my  advisor.  Dr.  Charles  Bridgman  of  the  AFIT 
Physics  Department.  I  hope  this  report  demonstrates  that 
his  encouragement  and  faith  in  me  has  been  rewarded. 

I  would  also  like  to  thank  Mrs.  Marge  Ilockemeier,  wife 
of  my  classmate  John,  for  typing  much  of  my  draft.  Finally, 

I  would  like  to  express  my  appreciation  to  my  wife  Bonnie 
for  putting  up  with  me  during  this  hectic  period. 
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Abstract 

This  report  presents  a  code,  developed  by  the  author, 
that  predicts  neutron  and  gamma  fluences  (or  doses J,  including 
neutron  induced  gammas,  in  an  exponential  atmosphere  for 
burst  altitudes  between  7  and  100  kilometers.  The  problem 
was  solved  by  using  the  diffusion  equation,  with  separation 
of  virgin  and  scattered  particles,  in  a  .lonorthogonal  coordi¬ 
nate  system.  The  diffusion  equation  in  this  coordinate 
system  was  approximated  by  a  nine  point  difference  equation 
and  the  resulting  matrix  equation  was  solved  by  use  of  a 
block  tri-diagonal  algorithm.  The  resulting  computer  code, 
in  FORTRAN  Extended,  was  written  to  calculate  the  surviva¬ 
bility  of  up  to  100  aircraft  or  space  vehicles  in  addition 
to  the  calculation  of  fluences  (doses).  The  code  requires 
20  minutes  of  central  processor  time  and  one  hour  of 
input/output  time  on  the  Wright-Patterson  Air  Force  Base 
(!DC  6600  computer.  The  results  of  this  code  for  a  burst  at 
25  km  are  compared  to  those  obtained  from  a  constant  density 
atmospheric  model  and  from  charts  based  on  SMAUG  which 
employs  mass  integral  scaling.  Significant  differences  are 
noted  for  the  case  where  the  burst  and  receiver  are  at  the 
same  altitude,  which  casts  some  doubt  on  the  validity  of 
mass  integral  scaling  at  this  altitude. 
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A  CODE  FOR  AIRCRAFT  SURVIVABILITY 
ANALYSIS  -  GAMMA  AND  NEUTRON  EFFECTS 

I .  Introduction 

Numerous  attempts  have  been  made  to  determine  the 
radiation  field,  neutrons  and  gammas,  from  nuclear  weapons. 
Straker,  in  a  report  titled  Status  of  Neutron  Transport 
in  The  Atmosphere  (Ref  5)  has  listed  these  attempts.  He 
lists  a  total  of  35  attempts  of  which  21  were  based  on  a 
constant  density,  (no  variations  with  al:itude,  no  ground, 
no  clouds)  infinite  air  model.  Five  attempts  were  based  on 
an  exponential  variation  in  air  density  and  nine  were  con¬ 
cerned  only  with  the  air/ground  interface.  The  method  of 
solution  for  26  of  the  35  attempts  was  Monte  Carlo  computer 
codes.  All  of  the  exponential  air  models  were  solved  by 
Monte  Carlo  codes.  Monte  Carlo  codes  normally  require  hours 
of  computer  run  time  and  are  therefore  relatively  expensive 
calculations . 

Recently,  two  high  speed  computer  codes  have  been 
developed  to  estimate  the  neutron  and  gamma  radiation  dose 
in  the  vicinity  of  an  atmospheric  nuclear  detonation  by 
interpolation  of  a  library  of  precalculated  results.  These 
precalculated  results  include  some  of  those  listed  by 
Straker  and  others  by  the  code  authors.  Both  codes  were 
presented  at  the  Radiation  Transport  in  Air  seminar  held 
at  Oak  Ridge  National  Laboratory  nn  15-17  November  1971. 

One  code,  ATR  (Ref  10),  was  developed  by  Science  Applications, 
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I:n.c.  The  final  documentation  of  this  code,  however,  has 
not  '/it  been  made  available  to  the  Radiation  Shielding 
Eh  font  ation  Center  (RS1C)  of  Oak  Ridge  National  Laboratory. 
The  ether  code,  SMAUG  (Ref  4),  was  devsloped  by  Murphy  of 
the  Air  Force  Weapons  Laboratory  (AFWL) ,  Kirtland  Air  Force 
3a.se,  New  Mexico.  SMAUG  calculates  neutron  and  secondary 
;amma  fluences  through  a  series  of  equations  that  were  curve 
fitted  to  the  data  of  Straker  and  Gritzner  (Ref  6).  The 
calculation  of  primary  gamma  fluences  is  done  through  a 
series  of  equations  that  were  curve  fitted  to  the  data  of 
3igoni  (Ref  4:4).  Both  data  sources  were  generated  using 
a  homogeneous  constant  -  density,  infinite  air  model. 

Both  SMAUG  and  ATR  require  only  seconds  of  computer 
time  (on  a  CDC  6600)  and  thus  both  offer  inexpensive  air 
transport  calculations.  However  they  both  rely  on  a  library 
of  homogeneous  air  calculations  which  are  modified  by  mass 
integral  scaling  to  account  for  the  true  exponential  varia¬ 
tion  of  the  air.  This  approach  is  probably  valid  only  to 
burst  altitudes  of  20-25  km.  The  code  ATR  accounts  for  the 
air-ground  interface  by  the  first-last  collision  approxima¬ 
tion  (Ref  9) . 

Similar  inexpensive  calculations  are  needed  for  burst 
altitudes  of  2S  km  and  above,  particularly  for  ABM  war  gaming 
and  reentry  vehicle  survivability  calculations.  Such  a 
technique  is  reported  here.  The  author  has  sacrificed  the 
accuracy  of  Monte  Carlo  and  higher  order  Boltzmann  equation 
solutions  as  used  in  the  3S  previous  attempts  reported  by 
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Striixcr  in  favor  of  the  Diffusion  (P-1)  approximation  of  the 
Boltzmann  Equation.  However  some  of  the  usual  weaknesses  of 
this  approximation  (diffusion)  are  avoided  by  a  mathematical 
separation  of  scattered  and  unscattered  radiation.  The 
latter  is  calculated  rigorously  and  the  diffusion  approxima¬ 
tion  is  used  only  for  the  scattered  radiation.  The  exponen¬ 
tial  variation  of  the  atmosphere  is  treated  with  a  unique 
coordinate  system  developed  here  by  the  author  which  also 
would  permit  the  inclusion  of  layered  clouds.  This  last 
feature  is  not  included  in  the  code  and  sample  calculations 
presented  here. 

The  resulting  code  is  designed  to  operate  as  a  sub¬ 
routine  of  a  nuclear  survivability  code  being  simultaneously 
developed  at  APIT  by  DeRaad  (Ref  2).  As  such,  the  input 
includes  the  spatial  position  of  up  to  100  target  vehicles 
and  their  vulnerabilities  to  neutrons  and  gammas.  The  code 
compares  the  radiation  field  at  each  vehicle  location  to 
the  vulnerability  level  for  a  survivability  determination. 
Additionally  the  code  will  output  iso-fluence  lines  for 
neutrons  and  gammas  including  neutron-induced  gammas. 

While  the  code  does  treat  the  exponential  air  exactly, 
the  goal  of  an  inexpensive  calculation  was  not  realized. 

A  calculation  requires  20  minutes  of  central  processor  (CP) 
time  and  60  minutes  of  input-output  time  on  the  Wright- 
Patterson  Air  Force  Base  CDC  6600  computer.  While  this  is 
less  than  a  Monte  Carlo  calculation  it  is  still  excessive 
for  repeated  runs.  (Perhaps  $400  per  run).  However  the 
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remits  do  offer  an  interesting  test  of  the  nass  integral 
scaling  approximation  against  a  true  exponential  atmosphere 
a.T.-i  provide  soao  insight  to  the  radiation  fields  from  high 
3.1  titjde  bursts. 

‘[he  mathematical  development  of  the  code  is  presented 
l  in.  Ch4pter  II.  The  use  of  the  code  is  illustrated  with  a 
!  s.npl<  problem  in  Chapter  III.  This  chapter  also  serves  as 
-i  complete  users  guide  to  the  code.  Chapter  IV  discusses 
■the  results  obtained  from  this  sample  problem.  Conclusions 
are  drawn  on  the  validity  of  these  results  and  the  success- 
fulness  of  the  code  in  Chapter  V.  Recommendations  are  also 
given  in  Chapter  V.  The  code  is  presented  in  Appendix  B. 
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XI.  Mathematical  Development 

The  development  of  this  computer  code  required  a 
numerical  approximation  of  the  diffusion  equation.  The 
numerical  approximation  depends  upon  the  atmospheric  model 
selected  and  the  coordinate  system  used.  Therefore,  the 
atmospheric  model  selected  is  discussed  first.  Next,  the 
diffusion  equation  is  presented  with  a  discussion  of  the 
coordinate  system  selected.  The  actual  numerical  approxima- 
tion  used- is  then  developed.  This  is  followed  by  a  discussion 
of  the  meshing,  source  terms,  boundary  conditions,  and  cross 
sections  used.  Finally,  the  actual  method  of  solution  is 
presented. 


The  Atmospheric  Model 

The  atmospheric  model  chosen  is  based  on  four  assump- 
tions.  The  first  assumption  is  that  the  composition  of  the 
atmosphere  is  21%  oxygen  and  79%  nitrogen.  The  second 
assumption  is  that  the  total  particle  density  varies  exponen¬ 
tially  according  to  the  relation 


(1) 


where 

p  »  particle  density  at  altitude  z  (particles/cm3) 
Pq  *  particle  density  at  sea  level  (particle/cm3) 

Z  ■  altitude  (km) 

H  •  atmospheric  scale  height  (km) 
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In  order  to  evaluate  the  constants  and  H,  this  equation 
was  changed  to  linear  form  by  taking  *he  natural  logarithm 
of  both  sides  to  get 

Anp  =  -Z/H  ♦  £np0  (2) 

t 

A  curve  fit  was  established  using  the  values  of  particle 
densities  from  the  U.  S.  Standard  Atmosphere,  1962  (Ref  7: 
2-19)  for  altitudes  of  0  to  100  km.  The  scale  height  H  was 

determined  to  be  7.0239  km  and  p^  was  determined  to  be 

19  3 

3.066  x  10  particles/cm  . 

The  third  and  fourth  assumptions  are  that  no  atmosphere 
exists  above  100  km,  and  that  the  earth  is  flat  throughout 
the  region  of  interest. 

The  Diffusion  Equation 

The  multigroup,  time  independent  diffusion  equation  for 
any  energy  group  g  is 

V-Dg(r)VFg(r)  -  Zg(r)Fg(r)  +  Sg(7)  =  0  (.3) 

K 

where 

(7  — 

D6(r)=  group  diffusion  coefficient  (cm)  as  a 
function  of  spacial  position  r 

Q  — 

F  (r )=  group  neutron  (or  gamma)  fluence  (particles/ 

2  — 
cm  )  as  a  function  of  spacial  position  r 

ff 

Er  =  group  macroscopic  removal  cross-section 
K 

(cm  *)  as  a  function  of  spacial  position  r 
Sb  =  group  neutron  (or  gamma)  source  (particles/ 
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The  Diffusion  Equation  in  a  Cylindrical  Coordinate 
System.  If  the  diffusion  equation  is  expressed  in  a 
symmetric  cylindrical  (r,Z)  coordinate  system,  Eq  (3)  can 
be  written  as  r.  partial  differential  equation  in  r  and  Z 
(Ref  8:6).  The  equation  is 

Dg(z)  32Fg(r,Z)  +  Dg(Z)  3Fg(r,Z)  +  3Dg(Z)  3Fg(r,Z) 

3v2  r  3r  3Z  3Z 

+  Dg(z)  1  fS-^-,-2)  -  Zg(Z)Fg(r,Z)  =  -Sg(r,Z)  (4) 

3  Z 

When  this  equation  is  approximated  by  a  difference  equation 
the  spacing  between  adjacent  radial  points  can  vary;  but, 
this  radial  spacing,  oi.re  fixed,  can  not  vary  with  altitude 
In  order  that  an  accurate  determination  of  fluence  can  be 
made,  the  radial  spacing  should  ha/e  at  least  two  or  three 
mesh  points  per  neutron  (or  gamma)  mean  free  path.  The 
horizontal  mean  free  path,  however,  is  dependent  on  the 
altitude.  The  exact  dependence  can  be  derived  from  the 
following : 

Let  us  combine 

Ig(Z)=  p(Z)ug  (5) 


where 


a  _  ] 

l  (Z)=  group  macroscopic  cross  section  (cm  ) 


at  altitude  Z 
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3 

p(Z)  =  particle  density  (particles/cm  )  at 
altitude  Z 

o  .  .  ,  -  2 . 

0°  =  group  microscopic  cross  section  (cm  ) 

with  Eq  (1)  to  obtain 

<7  o  -7  /  Jt 

Z'-'(Z)  =  (Z  =  0) e  '  (6) 


where 


Zg(Z=0) 


group  macroscopic  cross  section  (cm  ) 
at  sea  level 


The  horizontal  mean  free  path  is  dependent  on  the  macroscopic 
total  cross  section  and  is  given  by  the  relation 

Xg  =  1/Zg  (7) 

where 

Ag  =  group  mean  free  path  (cm) 

Eg  =  group  macroscopic  total  cross  section  (cm  *) 

The  dependence  on  altitude  Z  can  be  shown  by  combining 
Eqs  (7)  and  (6)  to  get 


Ag 


eZ/» 

Eg(Z=0) 


(8) 


Equation  (8)  clearly  illustrates  that  the  horizontal  mean 
free  path  increases  exponentially  with  altitude.  The 
required  radial  spacing  is  therefore  controlled  by  the 
lowest  altitude  of  interest.  This  requirement,  however, 
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neans  that  the  radial  spacing  at  the  highest  altitude  of 
interest  will  be  a  very  snail  fraction  of  a  mean  free  path. 

Since  a  difference  equation  is  written  for  each  radial 
point  on  each  altitude  line,  the  number  of  points  should  be 
a  minimum.  However,  at  the  higher  altitudes,  more  radial 
points  exist  than  required  for  a  sufficiently  accurate 
solution.  This  observation  implies  that  another  coordinate 
system  may  be  advantageous.  Since  the  solution  requires 
about  three  points  per  mean  free  path,  one  coordinate 
logically  should  be  mean  free  path. 

A  Non-Ortho gon a 1  Coordinate  System.  The  author  there¬ 
fore  proposes  the  following  coordinate  transforms: 


1  „  1  eA  2 

y  =  Hx  e  cos  x 


2  „  1  eA  .  2 

y  =  Hx  e  sxn  x 


(10) 


3 

y  = 


(ii) 


where 

y*  =  the  x  coordinate  in  the  Cartesian  coordinate 
system 

2 

y  =  the  y  coordinate  in  the  Cartesian  coordinate 
system 

y3  =  the  z  coordinate  in  the  Cartesian  coordinate 


system 


GNE/PH/72-8 


H  =  the  scale  height  of  the  atmosphere. 

x*  =  the  first  new  coordinate 
2 

x  =  the  second  new  coordinate 
3 

x  =  the  third  new  coordinate 


x  ,  x  ,  and  x  are  defined  to  be 

x1  -  N 

HE|(Z=0) 


x3  =  £n (Z/H) 


(12) 


(13) 


(14) 


where 

N  =  the  radial  distance  in  mean  free  paths  from 
the  x3  axis 

3 

0  =  the  transverse  angle  about  the  x  axis 


The  diffusion  equation  now  becomes 


„  n  r  1  2  5.nr  <  12  3..  _  ,  1  2  3.„,  12  3.. 

V-D(x  , x  , x  )VF(x  ,x  ,x  )  -  Er(x  ,x  ,x  )F(x.  ,x  ,x  ) 


12  3 

+  S(x  ,x  ,x  )  =  0 


(15) 


in  this  new  coordinate  system.  In  tensor  notation,  the  first 
term  of  Eq  (15)  can  be  written 


V-DVF.-i  l  l 

✓g  K=1  3x  l  P=1  3x‘ 


(16) 
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where 

KP 

g  are  the  contravariant  metric  tensors 
g  is  the  determinant  formed  by  gKp>  the  metric 
tensors 

Therefore,  b^th  the  metric  and  the  contravariant  metric 
tensors  need  to  oe  determined  to  evaluate  Eq  (16). 

The  metric  tensor  is  defined  as 


3  3  .  c  „  B 

V  r  3y  3y  r 

KP  L.  L.  .  K  .  P  °aB 
a=l  6=1  3x  3x 


L 


(17) 


where 


6  _  =  the  Kronecker  delta 
aB 


When  Eq  (17)  is  solved,  the  following  nine  metric  tensor 
components  are  obtained: 


gn  =  H  e 


2  2e 


(18) 


g 


12 


0 


(19 


3  x“ 
„2  1  x°  2e* 
g13  =  H  x  e  e 


(20) 


g2i  =  o 


(21) 


..2,  1.2  2e' 
g22  =  H  (x  )  e 


(22) 


-w- 
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g23  “  ° 


(23) 


3  x' 
„2  1  x  2e 
g31  =  H  x  e  e 


(24) 


g32  ~  0 


(25) 


u2  2x“  ri  ,  ,  1.2  2e"  . 
g33  =  H  e  [1  +  (x  )  e  ]  (2b) 


Therefore*  g  is 


g  = 


11 

g12 

g13 

21 

g22 

g23 

u6  (  12  2xJ  4eX 
=  H  (x  )  e  e 

(27) 

31 

g32 

g33 

KP  . 


The  contravariant  metric  tensor  g  is  defined  as 

.KP 


g 


KP  G 


g 


where 


(28) 


KP 

G  =  the  cofactor  of  the  element  gj,p  in  the 


determinant  of  Eq  (27) 

When  Eq  (28)  is  solved,  the  following  nine  contravariant 
metric  components  are  obtained: 


i  ,  1.2  2e' 

,11  1  +  (x  )  e 


u2  2e‘ 
H  e 


(29) 
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12  n 
g  =0 


13  -  x 


2  x' 
H  e 


g21  -  0 


22 


u2  f  L2  7e 
H  (x  )  is 


23 

g  =  0 


31 


-x 


u2  x* 
H  e 


32  n 
g  =0 


33 


u2  2x' 
H  e 


(30) 


(31) 


(32) 


(33) 


(34) 


(35) 


(35) 


(37) 


If  we  assume  fluence  is  symmetric  with  respect  to  x  and 

KP 

expand  Eq  (15),  eliminating  the  zero  values  of  g  ,  Eq  (16) 
becomes 


V • DVF  =  — 

/? 


_9 

3x 


r  n  /  n  3F 


g 


1. 3  9F_ 
9x 


9  r  _  /  31  9F  33  9F 

♦  — v  /g  D  f  g  — r  +  g  - t 

9x  \  9x  9x 


(38) 
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In  order  to  evaluate  this  expression,  the  functional  depend¬ 
ence  of  D,  the  diffusion  coefficient,  must  be  determined. 

The  diffusion  constant  is  defined  as 


D 


3^(Z) 


(39) 


where 


Z^.^(Z)  *  the  macroscopic  transport  cross  section 
(cm  * )  at  altitude  Z 


The  transport  cross  section  is  defined  as  £_  =  E  -  y£ 

IK  X.  o 

where  y  is  the  average  cosine  of  the  angle  of  anisotropic 
scatter  and  Eg  is  the  macroscopic  scatter  cross  section. 
Therefore,  combining  Eq  (6)  and  Eq  (39),  D  becomes 


When  Eq  (11),  the  Z  coordinate  transform  equation,  is 
used  on  Eq  (40),  D  becomes  a  function  only  of  x  and  is 
given  by 

3 

x 

D  =  D0ee  (41) 

If  the  expressions  for  D,  g,  g11,  g13,  g31,  and  g33 
are  substituted  into  Eq  (38)  and  the  indicated  partial 
derivatives  evaluated,  Eq  (38)  becomes 
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1  3F  ,  eX  -  1  3F 

3  .  1  +  ,3  a  3 

,  -  x  3x  2x  ox 

1  2e  e 

x  e 


32F 


3xX3x 


3 


.  A  ,  1,2  2e' 
1  +  (x  )  e 


2e 


32F 


ar  1,2 

3  (x  ) 


2x‘ 


32F 


»C*V 


-  5^  VZ=0>F  -  - 


H 


Doe 


(44) 


The  Difference  Form  of  the  Diffusion  Equation 

In  order  to  express  Eq  (44)  in  difference  form,  the 
dependent  variables  Fg  and  Sg  at  any  point  (x^,  x^ )  are 
defined  as 

Fg(xf,  x3)  =  F?  .  (45) 

i  3  i,3 

S8(x? ,  x3)  =  S?  .  (46) 

i  J  1,3 

The  partial  derivatives  of  Eq  (44)  can  be  expressed  in 

difference  form  by  the  use  of  a  central  difference  operator. 
The  partial  derivatives  are  therefore 


3F?  .  F2  .  -  F2  . 

_ LdL  =  _ izLi  l 

t  1 

3x'  2 Ax 


(47) 
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—  + 


1.2  2eJ 
(Ax  )  e 


Ax 


F? 


on  1  1  2~ 

2Ax  x  e 


-3  .  i  +  l>3 


♦  - ^ - ,  F? 


i  \ 


1  7  3  i-l.j-l  -  O  O  3  ,3 

1.  3  x  l  ,,  o.2  2x  ..  3  x 

2Ax  Ax-  e  V  (Ax  )  e  2Ax  e 


3  2x' 
2Ax  e 


F?  . 


H2S? 


i  ,j  +  l 


F. 


1 


2Ax  Ax  e 


3  x3  i+1»j 


+  1 


i  >3 


D?ee 


(52) 


The  fluence  coefficients  of  Eq  (52)  are  defined  as 


G.  .  =  - - - : 

1,3  ..  1.  3  x' 

2Ax  Ax  e 


(53) 


A.  = 


1 


3 


1  1 
- -  +  - 


1  (Ax3)2e2x  2Ax3eX~  2Ax3e2x 


(54) 


B.  . 


L2  2e' 
(Ax  )  e 


1  1  2e' 
2Ax  x  e 


(55) 


C.  - 
i,3 


or 

2  (x  ) 

+  —  -V-  -a  + 


, .  1.2  2e 
(Ax  )  e 


X"  (Ax3)2e2x 


o  n2£p (Z=0) 

- - 7  +  - - - (56) 


D 


0 


B.  . 
i,3 


, .  1.2  2e‘ 
(Ax  )  e 


oa  1  i  2eX 

2Ax  x  e 


(57) 
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A.  =  - i - = 

3  rA  3,2  2x' 
(Ax  )  e 


(58) 


^  *  3  x* 
2Ax  e 


3  2x' 
2Ax  e 


When  Eqs  (53)  through  (58)  are  substituted  into  Eq  (52), 
the  difference  equation  finally  becomes 


-G.  .F? 


+  A.  F? 


+  G.  .F?  ,  .  +  B.  .F? 


i,j  1-1, j-l  3  1,3-1  1,3  l+l, 3-1  1,3  1-1,3 


-C.  . F?  .  +  B.  . F? 


+  G.  .Ff 


+  A.F? 


i,j  i ,j  i , j  i+l,j  i,j  i-l,j+1  3  i,3  +  3 


■G.  .  F?  .  , 
1,3  l+l, 3+l 


K2S?  . 
1.3 


Mesh,  Source ,  Boundary  Conditions ,  and  Cross  Sections 

In  order  to  completely  prepare  Eq  (59)  for  use  in  a 

computer  four  factors  must  be  defined.  First  the  mesh  area 

1  3 

and  mesh  interval  (that  is,  the  range  of  x  and  x  and  the 
1  3 

values  of  Ax  and  Ax  )  must  be  determined.  Second,  the 

CT 

source  term  S?  ^  must  be  defined.  Third,  boundary  condi¬ 
tions  around  the  edge  of  the  mesh  must  be  defined.  Fourth, 
the  cross  section  used  must  be  defined. 

The  Mesh  Area  and  Mesh  Interval .  The  x*  coordinate  is 
a  function  of  the  sea  level  mean  free  path  of  some  energy 
group  of  gammas  or  neutrons.  This  can  be  shown  by  combining 
equations  (12)  and  (7)  to  get 


1  NAg(Z=0) 
X  "  H 
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where 

A^(Z=0)  =  the  sea  level  horizontal  mean  free  * 
of  gammas  or  neutrons  of  an  energy 
group  g. 

X  3 

The  minimum  value  of  xA  is  zero  (the  x  axis)  - 
Rigorously,  the  maximum  value  of  x*  should  be  infinite 
because  the  neutrons  and  gammas  are  transported  to  great 
heights  or  radial  distances.  However,  because  a  finite 
number  of  mesh  points  must  be  defined  for  a  computer  solu¬ 
tion,  these  boundaries  must  be  defined  as  a  finite  (and 
large)  number  of  mean  free  paths,  N,  and  the  fluence  at 
these  boundaries  will  be  forced  to  zero.  In  order  to 
determine  the  effect  of  this  approximation,  the  author  wrot 
a  homogeneous  air,  constant  density,  one  energy  group  code 
in  which  the  boundary  was  first  taken  as  10  mean  free  paths 
and  then  as  25  mean  free  paths.  The  fluences  were  identical 
up  to  seven  mean  free  paths  from  the  burst.  Between  seven 
and  nine  mean  free  paths  only  the  third  significant  figure 
of  the  fluences  were  different.  Only  in  the  tenth  mean 
free  path  order  of  magnitude  differences  appeared.  Therefore 
the  maximum  value  of  x*  was  chosen  to  correspond  to  ten  mean 
free  paths. 

The  value  of  A®,  the  mean  free  path,  could  be  selected 
from  any  energy  group.  However,  in  order  to  obtain  reasonabl 
accuracy  the  maximum  mean  free  path  was  chosen. 
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Ideally,  the  mesh  interval  should  be  chosen  such  that 
the  shortest  aean  free  path  has  at  least  three  aesh  points . 
This  requires,  in  the  40  group  structure  used,  about  600 
total  mesh  points  over  the  x*  axis.  This  nunber  of  points 
was  inpractical  to  use  on  the  conputer  available.  Therefore. 
40  points  were  selected  for  the  sample  run  of  Chapter  III. 
(However,  the  program  was  constructed  so  that  up  to  60  points 
could  be  used.)  The  description  of  the  MESH  subroutine  in 
Appendix  B  describes  how  this  option  may  be  used. 

The  determination  of  the  range  of  x"5  was  also  based  on 
a  maximum  range  of  ten  mean  free  paths  from  the  burst.  In 
this  case  the  maximum  value  of  xJ  would  be  equivalent  to 
ten  mean  paths  up  from  the  burst  and  the  minimum  value  would 
be  equivalent  to  ten  mean  free  paths  down  from  the  burst. 

The  energy  group  with  the  longest  mean  free  path  was  also 

3 

chosen  to  evaluate  the  range  of  x  .  However,  a  simple 

3 

relationship  between  x  and  mean  free  path  such  as  Eq  (60) 
does  not  exist  since  the  air  density  is  changing  over  the 
path  traveled  by  the  neutron  or  gamma.  A  relationship  can 
be  derived  between  altitude  Z  and  mean  free  path  based  on 
the  assumption  of  the  exponential  variation  in  air  density 
with  height. 

This  relationship  is  based  upon  the  definition  of 
macroscopic  cross  section  given  in  Eq  (5): 

Et(Z)  =  P(Z)aT  (5) 
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If  this  equation  is  integrated  over  the  average  mean  free 
path  at  any  altitude  2  to  an  altitude  Z  +  X 


Z+X 


Z+X 


fz  Zt-'  *  f  «!**’ 


(61) 


the  following  series  of  relationships  can  be  obtained. 

On  the  left  hand  side  of  Eq  (60)  we  extract  the  average 
vali.e  of  E^  over  one  mean  free  path  (  <E^>  )  and  on  the  right 
hand  side  replace  p  by  its  definition  in  Eq  (1)  to  get 


Z+X 


<£t>  sz  «•  ■  st 


z+x 


-Z*/H._t 
aTP0e  <xl 


(62) 


But  the  integral  on  the  left  hand  side  is  X  and  the  expression 
pQaT  on  the  right  hand  side  is  Et(Z=0),  a  constant.  Therefore 


Z+X 


<£t>  X  =  Et(Z=0)  J  e 


-Z'/H 


dZ' 


(63) 


Carrying  out  the  integration  on  the  right  side  and  using 
the  definition  of  mean  free  path  as  the  reciprocal  of  macro¬ 
scopic  total  cross  section,  Eq  (63)  becomes 


.  _  -Z/H..  -X/H 

1  =  Et(Z=0)e  il  -  e 


But,  Et(Z-0)e  is  (Z) .  Therefore 


] 


(64) 


1  =  EtU  -  e'X/H] 


(65) 


Solving  Eq  (65)  for  \ 
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(66) 


Therefore,  this  is  the  relationship  between  a  and  Z  (since 
is  a  function  of  Z)  for  the  average  rean  free  path  up 
fron  any  altitude  3. 

In  a  sinilar  nanner,  the  relationship  between  X  and  Z 
for  the  average  nean  free  path  down  fron  any  altitude  Z  can 
be  found  to  be 


(67) 


Equations  (66)  and  (67)  were  used  to  find  the  upper  and 
lower  values  of  altitude  equivalent  to  ten  mean  free  paths 
in  either  direction.  The  minimum  and  maximum  values  of  x3 
were  then  obtained  from  its  definition  of 

=  in (Z/H)  (14) 

3 

The  minimum  value  of  x  ,  however,  has  an  additional 
constraint.  Z  can  not  be  equal  to  or  less  than  zero.  This 
constraint  in  turn  causes  a  lower  limit  of  seven  kilometers 
for  the  burst  height  in  order  to  keep  a  ten  mean  free  path 
lower  boundary  without  intersecting  the  flat  earth.  The 
adoption  of  an  upper  limit  of  100  km  for  the  model  atmosphere 
places  a  restriction  upon  the  upper  value  of  x3  to  100  km. 

As  in  the  case  for  the  x*  mesh  interval,  the  mesh 
interval  should  be  selected  such  that  the  shortest  mean  free 
path  has  at  least  three  mesh  points.  This  requires,  in  the 
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40  group  structure  used,  about  2200  total  aesh  points  over 

3 

the  x  axis.  Again  to  limit  the  number  of  mesh  points  in 
„  3 

tne  x  direction,  70  points  were  selected  for  the  sample 
run  of  Chapter  III.  (The  program,  however,  was  constructed 
to  allow  up  to  120  points  in  the  x-5  direction.)  The 
description  of  the  MESH  subroutine  in  Appendix  B  describes 
how  this  option  say  be  used. 

Definition  of  the  Source  Tern.  The  diffusion  equation 
inplicitly  assunes  a  linear  variation  in  the  cosine  of  the 
directional  angle  of  the  diffusing  particles.  This  follows 
because  the  diffusion  equation  results  frora  a  Legendre 
expansion  in  direction  of  the  sore  rigorous  Boltzmann 
equation  which  retains  only  the  first  two  terras.  .Vear  the 
source  the  virgin  particles  make  up  the  largest  fraction  of 
the  total  fluence  and  are  nearly  monodirectional  in  the  out¬ 
bound  direction.  Thus  diffusion  theory  is  weakest  close  to 
the  source  which  is  the  region  of  highest  interest.  This 
dilemma  can  be  avoided  by  defining  the  fluence  at  any  point 
to  have  two  components. 

The  first  component  at  any  point  consists  of  the  virgin 
(unscattered)  particles  from  the  burst.  The  second  component 
consists  of  the  scattered  particles.  This  second  component 
at  any  point  therefore  contains  those  particles  scattered 
down  from  a  higher  energy  group  or  from  virgin  scatter  which 
entered  the  energy  group  at  some  other  point  and  have  not 
yet  downscattered  to  a  lower  energy  group. 
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The  virgin  particles  are  solved  rigorously  while  the 
diffusion  equation  is  used  to  describe  only  the  scattered 
particles.  The  source  term  in  this  diffusion  equation 
consists  of  those  particles  scattering  out  of  the  virgin 
groups  at  that  point  plus  those  particles  downscattering 
from  higher  energy  scattered  groups.  The  source  term  can 
therefore  be  expressed  as 


s' 

i  Ti 

g*=i  j 


F 


g' 

i»  5 


(68) 


where 


Fg’ 

v. 

i.-l 


=  the  source  term  for  group  g  at  spacial 

3 

point  i,j  (particles/ciu  ) 

=  the  macroscopic  scatter  cross  section  for 

virgin  group  g'  (cm  1) 

=  the  virgin  fluence  in  group  g'  at  spacial 

2 

location  i,j  (particles/cra  ) 

=  the  macroscopic  scatter  cross  section  for 

scattered  group  g'  (cm  *) 

=  the  scattered  fluence  in  group  g'  at 

2 

spacial  location  i,j  (particles/cm  ) 


However,  in  order  to  evaluate  Eq  (68),  the  virgin 
fluence  for  any  point  i,j  must  be  determined.  The  virgin 
fluence  at  any  point  is 


si  e 


-at->  r 


v .  - 


(69) 


4tiR‘ 
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where 

% 

g 

S*  =  the  total  number  of  particles  in  group  g 
emitted  from  the  nuclear  weapon. 

R  =  the  distance  from  the  burst  point  to  the 
point  i,j 

=  the  average  microscopic  total  cross  section 
of  the  particle  over  the  distance 

The  value  of  can  found  by  integrating  the 

point  value  of  (which  is  dependent,  upon  the  altitude  Z) 
over  the  path  length  R.  So 

R 

J  lt  dR*  =  <Zt>  R  (70) 

However,  since  1^.  is  a  function  of  altitude  the  left 
side  should  be  expressed  in  terms  of  Z.  Figure  1  illustrates 
the  geometry. 


HOB 


AZ 


Z 


Fig.  1.  Sample  Geometry  Relating  Variables  Needed  in 
Calculation  of  Average  Cross  Section.  (Note:  HOB  is 
Height  of  Burst.) 
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From  Fig.  1  the  following  relationship  is  obvious: 

\ 

dR  =  esc  6dZ  (71) 

Therefore,  from  Eqs  (71),  (70),  and  (6) 

2 

<E  >  R  =  Z  (Z  =  0)  esc  0  f  e'Z'/H  dZ  •  (72) 

Z  Z  *'H0B 


Once  this  integration  is  performed,  Eq  (72)  becomes 

-  <E  >  R  =  HZ  (Z=0)  esc  0[e-Z/H  -  e'K0B/H]  (73) 

»*  t 

and  also 

CSC  9  =  Bj  (74) 

Therefore,  from  Eqs  (69),  (73),  and  (74), the  virgin 
fluence  at  any  point  is  given  by 


Sq  HZt(Z=0) 
2  6 

4ttR 


R 

AZ 


-HOB/H. 
e  ] 


(75) 


Using  Eqs  (75),  (68),  and  (59),  a  difference  equation 

can  be  written  for  any  interior  point  i,j  in  the  meshed  area 

13  1 

defined  by  the  ranges  of  x  and  x  and  by  the  intervals  Ax 

3 

and  Ax  .  However,  since  Eq  (59)  is  a  nine-point  difference 
equation  involving  fluences  at  points  such  as  i-1,  j-1, 
boundary  conditions  must  be  defined. 

Boundary  Conditions.  The  first  boundary  condition  that 
must  be  defined  is  for  the  points  on  the  x^  axis,  that  is, 
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when  x*  is  zero.  Two  problems  appear  if  the  difference 

equation  (59)  is  used.  This  first  problem  is  that  the 

coefficients  B.  .  and  B.  .  contain  the  term  1/x*  which  is 
i,3  i,3 

infinite  at  x  equal  to  zero.  The  second  problem  is  that  at 

x*  equal  to  zero,  i  is  equal  to  one,  and  in  the  nine  point 

difference  equation  three  i  -  1  (i  =  0)  terms  appear. 

The  second  problem  is  easily  treated  since,  in  the 

1  2  3 

development  of  the  diffusion  equation  in  the  x",  x  ,  x 

coordinate  system,  symmetry  of  fluence  with  respect  to  the 
2 

x  coordinate  (the  transverse  angle)  was  assumed .  Therefore 


F 


0,3 


(76) 


•The  first  problem  requires  a  special  development  of  the 

difference  equation.  This  requirement  exists  since  the 

12  3 

diffusion  equation  in  the  x  ,  x  ,  x  coordinate  system, 

Eq  (44) ,  contained  the  term 


1 


x 


1 

e 


Equation  (76)  is  equivalent  to  stating 


The  term,  therefore 
ilowever,  the  use  of 


at  x*  =  0  becomes  0/0  which  is  undefined. 
L'Hospital’s  rule  on  this  term  results  in 
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x1-+Q 


Lid 

1  2ex  x'-i-O 
x  e 


o  I  •> 

3  F/3 (x  j 


2e 


(78) 


Therefore,  using  the  substitution  indicated  in  Eq  (78) 
into  Eq  (44),  forming  a  new  difference  equation,  and  using 
the  definition  of  Eq  (76),  the  following  special  form  of 
the  difference  equation  for  i  =  1  is  obtained: 


AF,  -  ,  -  CF,  .  +  BF  .  ♦  AF.  .  . 
1.3-1  1 ,3  2,3  l,j+l 


H2S.  . 

- ^4  (79) 


V 


where 


A  = 


1  ♦  _ L 


3.2  2x3  3  x3  ,.3  2x3 

(Ax  )  e  2Ax  e  2Ax  e 


(80) 


C  = 


—  + 


H‘Zr(Z=0) 


.  _  -  x^  3-2  2xj 

(Ax1) 2e2e  <Ax  )  e 


(81) 


B  = 


, .  1,2  2e 

(Ax  )  e 


(82) 


A  = 


, A  3-2  2x' 
(Ax  )  e 


3  x* 
2Ax  e 


3  2x' 
2Ax  e 


(83) 


The  second  boundary  that  must  be  defined  is  when  x3  is 
a  minimum  (the  lower  boundary  of  the  mesh).  This  lower 
boundary  was  chosen  far  enough  away  from  the  bu.rst  (10  mean 
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free  paths  of  the  group  with  the  longest  mean  free  path)  so 

that  the  value  of  fluence  is  at  least  seven  or  eight  orders 

of  magnitude  lower  than  the  fluences  existing  around  the 

burst.  Therefore,  along  this  boundary,  fluence  is  assumed 

to  be  zero.  Reviewing  Eq  (53)  through  (59)  this  condition 

implies  that  the  coefficients  of  the  F.  ,  .  ,,  F.  .  and 

r  l-l, j-1*  i,j-l 

F.  ,  .  .  terms  must  be  zero  when  j  =  1. 
i+l, 1-1 

The  third  boundary  that  must  be  defined  is  when  x  is 

a  maximum  (the  outer  boundary  of  the  mesh) .  This  boundary 

was  also  chosen  so  that  the  fluences  existing  at  this 

boundary  were  low  compared  to  those  around  the  burst. 

Therefore,  along  this  boundary,  fluences  are  assumed  to  be 

zero.  This  condition  implies  that  the  coefficients  in 

Eq  (59)  of  the  F.  .  .  F.  .  and  F.  ,  .  terms  must  be 

M  l  +  l,  i-l*  i  +  l,l’  i  +  l, 1  +  1 


zero  when  x  is  a  maximum. 

The  final  boundary  condition  that  must  be  defined  is 

3 

when  x  is  a  maximum  (the  upper  boundary  of  the  mesh) .  Two 

possible  conditions  can  exist.  First,  this  upper  boundary 

is  less  than  100  km  and  second,  this  upper  boundary  is 

100  km.  If  the  first  condition  exists,  the  upper  boundary 

is  ten  mean  free  paths,  of  the  group  with  the  longest  mean 

free  path,  away  from  the  burst.  Therefore,  by  the  same  logic 

of  the  preceding  two  paragraphs,  fluence  along  this  boundary 

is  zero  and  the  coefficients  of  the  F.  ,  .  ,,  F.  .  ,,  and 

i-l,  1  +  1  i,l  +  l 

F i + i  terms  of  Eq  (59)  must  be  zero. 

The  second  condition  occurs  when,  in  the  calculation 
for  the  upper  boundary,  altitudes  of  greater  than  100  km  are 
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obtained.  In  this  case,  the  upper  boundary  is  set  to  100  km. 
Therefore,  the  fluences  at  this  boundary  nay  not  be  very 
much  lower  than  about  the  burst  and  assuming  zero  fluences 
for  this  boundary  would  be  invalid.  However,  since  a  vacuum 
is  assumed  to  exist  above  100  km,  a  valid  boundary  condition 
is  that  the  return  current  is  zero.  Or  using  Fick's  Law  it 
can  be  stated  as 


t  F  D  3F 
J-  =  4  *  2  Ti  *  0 

d 


where 


J  =  return  current  (particles/ cm) 

D  =  diffusion  coefficient  (cm)  at  100  km 

2 

F  =  fluence  (particles/cm  ) 


This  equation  can  be  differenced  using  a  central 

difference  operator  and  solved  for  F.  .  ,  to  get 

i,3  +  l 


F  =  F 

i,j  +  l  i  ,  j - 1 


Ax  F .  . 
x  ,  1 


When  this  is  substituted  into  Eq  (59)  the  following 
special  form  of  the  difference  equation  is  obtained: 


HS.  . 

AF.  .  .  +  BF.  .  .  +  CF.  .  +  BF.  .  - - (86) 

i,j-l  1-1,3  i,3  i+l,3  x3 

V* 
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where 


A  = 


3.2  2x' 
(Ax  )  e 


(87) 


B  = 


rA  1 a2  2e 
(Ax  )  e 


1  1  2e' 
2Ax  x  e 


.  3  x* 

2D()Ax1eX  e6 


(88) 


C  = 


_  ♦  2^b.2  + 

.3  ,,..1,2 


H"Er(Z=0) 


(Ax 


.  ,  x~  (Ax  )  3,2  2x‘ 

1,2  2e  v  3  (Ax  )  e 

)  e 


+  _ I _ ,  ♦  _ 1 


1 


Z  r'1  x3  3  x3 

n  .  3  2xJ  eX  _n  2ex  __  2x  e 

DqAx  e  e  2^06  2DOe  6 


1,2  2e 
(Ax  )  e 


B  =  - 1 - ,  .  IXI  *  1 


1  1  2e' 
2Ax  x  e 


,  3  x' 

2D0Ax1eX  ee 


(89) 


(90) 


A  difference  equation  now  can  be  written  for  every  point 
in  the  mesh  with  all  ♦r.as  defined  except  for  the  cross 
sections . 
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Cross  Secti.ons .  The  cross  sections  selected  were 
obtained  from  the  Radiation  Shielding  Information  Center 
(RSIC)  of  Oak  Ridge  National  Laboratory  from  the  RSIC  Data 
Library  Collection  DLC-14  (Ref  3).  These  cross  sections  are 
multigroup  and  consist  of  22  neutron  groups  coupled  with 
18  gamma  groups.  Since  the  gamma  groups  were  coupled  with 
the  neutron  groups,  calculations  of  secondary  gamma  fluences 
(caused  by  neutron  scattering  in  air)  are  possible  in  addi¬ 
tion  to  the  calculation  of  primary  neutrons  and  gamma 
fluences.  These  cross  sections  were  selected  since  they 
were  used  for  Straker  and  Gritzner’s  constant  density, 
infinite  air  calculations  (Ref  6).  Therefore,  these  cross 
sections  were  indirectly  used  by  SMAUG  since  SMAUG  curve 
fits  Straker  and  Gritzner’s  data.  Use  of  these  cross 
sections  allows  comparison  between  the  code  presented  in 
this  report  and  SMAUG. 

These  cross  sections,  however,  are  for  an  air  density 

3 

of  l.ii  mg/cm  which  is  equivalent  to  about  1  km  in  the  U.S. 
Standard  Atmosphere,  or  to  about  2  km  in  the  curve  fitted 
exponential  atmosphere  used  in  this  report.  Since  the  entire 
mathematical  development  of  this  report  has  been  based  on  sea 
level  cross  sections,  these  cross  sections  were  extrapolated 
to  sea  level  using  Eq  (6).  In  addition,  the  cross  sections 
are  differential  cross  sections  expanded  in  six  terms  of  a 
Legendre  expansion.  Therefore,  removal  and  transport  cross 
sections  had  to  be  calculated.  The  energy  ranges  of  the 
neutron  groups  are  presented  in  Table  I  and  the  entrgy  ranges 
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Table  I 


Neutron  Group  Energy  Ranges 


Neutron  Group 
Number 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 


Enerjy  Range 
(MeV) 

12.2  to  15.0 
10.0  to  12.2 
8.19  to  10.0 
6.36  to  8.19 
4.97  to  6.36 
4.07  to  4.97 
3.01  to  4.07 
2.46  to  3.01 
2.35  to  2.46 
1.83  to  2.35 
1.11  to  1.83 
0.55  to  1.11 
0.111  to  0.55 
0.00335  to  0.111 
5.83E-4  to  3.35E-3 
1.01E-4  to  5.S3E-4 
2.90E-5  to  1.01E-4 
1.07E-5  to  2.90E-5 
3.06E-6  to  1.07E-5 
1.12E-6  to  3.06E-6 
4.14E-7  to  1.12E-6 
0.0  to  4.14E-7 


4 


(a) 


(a)  3.35E-3  reads  as  3.35  x  10 


(Ref  4:3) 
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of  the  gamma  groups  are  presented 
cross  seel  ions  used  are  presented 


in  Table  II.  The  actual 
in  Appendix  C. 


Method  of  Solution 

The  group  scattered  fluences  were  solved,  one  group  at 
a  time  starting  with  the  highest  energy  neutron  group 
(group  1) .  A  difference  equation  was  written  for  each  point 
in  the  mesh.  The  result,  for  any  group,  is  a  matrix  equation 
of  the  form 

A  F  =  S  (91) 


where 

A  =  the  matrix  consisting  of  the  terms  A,  B,  C,  A, 
B,  and  G 

F  =  the  array  of  unknown  scattered  fluences 
=  the  array  of  source  terms 


The  matrix  A  is  a  block  tri-diagonal  matrix.  Appendix  A 
defines  a  block  tri-diagonal  matrix  and  presents  the  algorithm 
used  to  solve  matrix  equation  (93). 

A  problem  arose  during  the  computer  solution  using  the 
coordinate  system  and  meshing  developed  in  this  chapter. 
Negative  fluences  occurred  in  the  upper  portions  of  the 

meshed  area.  This  was  a  result  of  the  radial  spacing  between 

12  3  1 

mesh  spaces.  In  the  x  ,  x  ,  x  coordinate  system.  Ax  is  a 

constant;  however,  the  actual  distance  between  mesh  points 
is  increasing  exponentially  with  increase  in  altitude.  There¬ 
fore,  in  the  upper  portion  of  the  mesh,  the  fluence  was 
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r«c  /ou  /-7-v  c 

urn./  *  11/  *  i-  v* 


Gang a  Group  Energy  Ranges 


Ganma  Group 
Xunber 


Energy  Range 
(MeV) 

8.0 

to 

t-4 

o 

• 

o 

6.5 

to 

8.0 

5.0 

to 

6.5 

4.0 

to 

3.0 

3.0 

to 

4.0 

2.5 

to 

3.0 

2.0 

to 

2.5 

1.66 

to 

2.0 

1.33 

to 

1.66 

1.0 

to 

1.33 

0.8 

to 

1.0 

0.6 

to 

0.8 

0.4 

to 

0.6 

* 

0.3 

to 

0.4 

0.2 

to 

0.3 

0.1 

to 

0.2 

0.05 

to 

0.1 

0.02 

to 

0.05 

(Ref  4:4) 


changing  rapidly  between  mesh  points  and  causing  the  negative 
fluences.  In  order  to  correct  this  fault,  x*  was  not  allowed 
to  exceed  the  equivalent  of  200  km  at  the  outer  boundary. 
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When  xA  exceeds  this  value,  a  different  coordinate  systea  is 

5 

used  starting  on  the  next  value  of  x  -  This  systea  was 


11  2 
y  =  x  cos  x 


(92) 


2  1.2 

y  =  x  sm  x 


(S3) 


3  „  x' 

y  =  He 


(94) 


1  2 

x  and  x  are  therefore  the  r  and  9  of  a  normal  cylindrical 
coordinate  system.  The  difference  equation  from  using  this 
coordinate  system  is 


AF.  .  .  +  BF.  ,  .  -  CF.  .  +  BF.  .  .  +  AF.  .  . 
i,3-l  1-1,3  i,l+l 


h2s.  . 
1  >  3 
D» 


(95) 


where 


A  = 


1 


1  1 
+  - 


3.2  2x3  ..  3  x3  ..3  2x3 

(Ax  )  e  2Ax  e  2Ax  e 


(96) 


B  = 


H 


H' 


1,2  0  1.  1 

(Ax  )  2x  Ax 


(97) 


C  = 


2H 


H‘Er(Z=0) 


(Ax1)2  3,2  2x 

v  '  (Ax  )  e 


(98) 


s  .  ,  .  «2 


. .  1,2  ,  1.  1 

(Ax  )  2x  Ax 


(99) 


GHE/PH/72-8 


A  = 


—  + 


3.2  2x' 
(Ax  )  e 


.  3  xJ 
2Ax  e 


_ ,  3  2x' 
2Ax  e 


(100) 


Since  this  difference  equation  is  used  only  in  the  upper  area 

of  the  mesh,  only  three  boundary  conditions  need  be  defined. 

They  are  for  the  x°  axis,  the  upper  boundary,  and  the  outer 

boundary.  The  boundary  conditions  are  the  same  as  defined 

3  — 

before.  Therefore,  for  the  x  axis,  A  and  A  remain  as 
defined  by  Eqs  (96)  and  (100).  B  becomes  zero,  and 


B 


4H 


* ,  1.2 
(Ax  ) 


(101) 


C  = 


4H 


H  IR(Z=0) 


,A  1.2 

(Ax  ) 


,A  3.2  2x' 
(Ax  )  e 


(102) 


For  the  upper  boundary,  the  same  two  possibilities  exist. 
The  boundary  can  be  less  than  1J0  km  or  it  can  be  100  km.  If 
the  upper  boundary  is  less  than  100  km.  A,  B,  C,  and  B  remain 
defined  as  in  Eqs  (96)  through  (99)  and  A  is  equal  to  zero. 

If  the  upper  boundary  is  100  km,  B  and  B  remain  as  defined  in 
Eqs  (97)  and  (99),  A  is  zero,  and 


A  = 


[K  3.2  2x~> 
(Ax  )  e 


(103) 
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C  = 


2H 


1>2 
(Ax  ) 


3.2  2x' 
(Ax  )  e 


H‘£r(Z=0) 

— 


.  3  2x3  ex 
DqAx  e  e 


3  x 
„  x  e 

2Doe-  e 


5  -  „•» 

on  2x^  e 

2Doe  e 


(104) 


For  the  outer  boundary,  A,  B,  C,  and  A  remain  as  defined 
as  in  Eqs  (96),  (97),  (98),  and  (100).  B  is  zero. 

The  mesh  interval  Ax°  remains  the  same  since  x  has  not 
been  changed.  The  new  Ax1  (actually  a  Ar)  is  the  radial 
equivalent  to  the  last  Ax*  of  the  old  coordinate  scheme. 

The  computer  code  was  written  for  the  Kright-Patterson 
Air  Force  Base  CDC  6600  computer  using  an  on-line  plotter. 

The  language  used  was  FORTRAN  Extended. 
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III.  A  Users  Guide  to  the  Code 

This  chapter  presents  a  guide  to  the  user  on  how  to  use 
the  code.  The  input  for  a  sample  problem  is  shown  as  an 
illustration.  Three  types  of  punched  cards  are  required; 
control  cards,  data  cards,  and  end  of  job  cards. 

Control  Cards 

Four  control  cards  are  required.  These  four  cards 
preceed  any  other  cards.  Entries  on  all  these  cards  start 
in  column  one.  These  cards  are  unique  to  the  Wright-Patterson 
Air  Force  Base  CDC  6600.  Other  computers  will  require  dif¬ 
ferent  control  cards. 

Card  No.  1: 

This  card  contains  the  users  initials,  time  to  run  the 

code,  core  size  in  octal,  problem  number  assigned  to  the  user 

by  the  computer  center,  name,  telephone  number,  and  class. 

EXAMPLE:  RDM, T3400, CM 120000.  T7 1067 5 .MCLAREN , 2533886 , 

GNE  72 

The  time  required  for  running  the  code  is  3400  seconds  and  a 
core  size  of  1200Q0o  is  required. 

O 

Card  No.  2: 

This  card  selects  the  FORTRAN  compiler. 

EXAMPLE:  FTN . 

Card  No.  3: 

This  card  loads  the  program  and  causes  it  to  go  into 
execution . 

EXAMPLE:  LGO . 
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Card  No.  4: 

This  card  signals  the  end  of  the  control  cards. 
EXAMPLE:  A  multipunch  7/8/9  in  column  one. 

Data  Cards 

The  computer  code  itself  follows  the  control  cards. 
After  the  code,  a  series  of  data  cards  are  required.  The 
data  is  input  in  either  an  I  format  or  an  E  format.  All 
numbers  are  right  justified  in  the  field  specified  by  the 


format  statement 

.  Examples  of  correct  and 

incorrect  use  are 

shown  below,  b 

represents  a  blank  space. 

All  entries  start 

in  column  one. 

-Format 

Correct 

Incorrect 

214 

bb22bbl8 

22bbl8bh 

bbb5- 100 

b5b~ 100b 

lx , 1PE1 1 . 4 

bbl . 2345E+04 
b- 1 . 2345E+04 

1 . 2345E  +  0  4bb 

Card  No.  1: 

This  card  contains  two  numbers  in  214  format.  The 
first  number  is  the  number  of  neutron  energy  groups  and  the 
second  number  is  the  number  of  gamma  groups.  This  card  is 
included  as  part  of  the  input  library  and  normally  requires  no 
preparation  on  the  part  of  the  user.  However,  if  the  user 
desires  to  use  a  different  cross  section  data  set  with  dif¬ 
ferent  number  of  groups,  this  card  must  be  replaced  with  a 
card  prepared  by  the  user  in  the  same  format. 

EXAMPLE:  bb22bbl8 
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Card  Wo.  2: 

This  is  actually  a  deck  of  cards  containing  the  multi- 
group  cross  sections.  The  format  for  each  card  is  lx,lP7E11.4. 
This  deck  is  included  as  part  of  the  input  library  and 
normally  requires  no  preparation  on  the  part  of  the  user. 

If  the  user  desires  to  use  different  cross  sections  ,  this 
deck  must  be  replaced  with  a  deck  written  in  the  same  format. 
The  cross  section  data  supplied  by  the  user  must  be  in  the 
following  order.  All  cross  section  data  for  the  highest 
energy  neutron  group  must  be  first.  This  is  followed  by  the 
data  for  the  rest  of  the  neutron  groups  in  the  order  of 
descending  energy.  Next,  the  highest  energy  gamma  group 
data  follows.  The  data  for  the  gamma  groups  must  also  be 
arranged  in  order  of  descending  energy.  Each  group  G  must 
contain  a  series  of  cross  sections  whose  number  must  be  three 
more  than  the  total  number  of  groups  (neutron  plus  gammas). 

In  any  group  G,  the  cross  sections  must  be  arranged  in  the 
following  order: 

1  is  transport 

2  is  removal 

3  is  total 

4  is  scatter  (G  to  G) 

5  is  scatter  (G-l  to  G) 

6  is  scatter  (G-2  to  G) 

7  is  scatter  (G-3  to  G) 

etc . 
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Every  care,  except  for  the  final  card  in  this  deck  nust 
have  sevea  combers  in  the  specified  format.  In  cases  where 
the  groop  to  groop  scatter  cross  section  is  zero  or  do  not 
exist  (example:  G-S  when  G  is  5)  a  zero  must  be  entered  on 
the  card.  The  cross  section  data  used  in  the  input  library 
is  listed  in  Appendix  C. 

Card  So-  — . 

This  is  also  a  deck  of  cards  and  contains  the  response 
functions  required  to  translate  the  sultigroup  neutron 
floence  to  silicon  dose  in  rads.  The  format  for  each  card  is 
iX,l?7£11.4.  This  deck  is  supplied  as  part  cf  the  input 
library.  If  the  user  desires  to  use  different  response  func¬ 
tions,  he  must  replace  this  deck.  The  number  of  values  in 
this  deck  must  correspond  to  the  number  of  neutron  groups 
specified  on  card  number  one.  The  response  functions  should 
be  arranged  in  order  cf  decreasing  energy  groups.  The 
response  functions  used  in  the  input  library  are  presented 
in  Table  III.  These  functions  were  obtained  from  the  SMAUG 
code. 

Card  No.  4: 

This  is  also  a  deck  of  cards  and  contains  the  response 
functions  required  to  translate  the  multigroup  gamma  fluence 
to  silicon  dose  in  rads.  The  format  is  the  same  as  for  card 
number  three.  This  deck  is  supplied  as  part  of  the  input 
library  and  if  the  user  wishes  to  use  different  response 
functions,  he  must  replace  this  deck.  The  number  of  values 
in  this  deck  must  correspond  to  the  number  of  gamma  groups 
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Table  III 


Neutron  Response  Functions 


Energy  Group 

For  Silicon  Dose 
in  rads 

For  Tissue  Dos 
in  rads 

1. 

9.3SE-10  (a) 

6 . 2545E-09 

2 

9.66E-10 

5 . 7744E-09 

3 

8.74E-10 

5 . 1525E-09 

4 

S.56E-10 

4.9958E-09 

5 

2.15E-10 

4 . 4861E-09 

6 

1.41E-10 

4 . 1 199E-09 

7 

1.02E-10 

4 . 0662E-09 

8 

8.4E-11 

3.4281E-09 

9 

7.7E-11 

3. 1542E-09 

10 

6.9E-11 

3 . 0504E-09 

11 

5.5E-11 

2 . 6187E-09 

12 

4.8E-11 

2.0131E-09 

13 

3.4E-11 

1 .2918E-09 

14 

o 

4 . 2594E- 10 

15 

0 

1.9581E-11 

16 

0 

3.6590E-12 

17 

0 

1 . 1533E- 12 

18 

0 

1 . 084  IE-  12 

19 

0 

1 . 5460E-  12 

20 

0 

2 . 6671E-  12 

21 

0 

4 . 3388E- 12 

22 

0 

8.2591E-12 

GNE/PH/72- 8 


specified  on  card  number  one.  The  response  functions  used 
in  the  input  library  are  presented  in  Table  IV.  These 
functions  were  obtained  from  the  SMAUG  code. 

Card  No.  S: 

This  card  is  the  same  as  card  number  three  except  the 
response  functions  are  for  tissue  dose  in  rads.  All  instruc¬ 
tions  are  the  same  as  for  card  number  three.  The  response 
functions  used  in  the  input  library  are  presented  in  Table 

III.  These  functions  were  obtained  from  the  SMAUG  code. 

Card  No.  6: 

This  card  is  the  same  as  card  number  four  except  the 
response  functions  are  for  tissue  dose  in  rads.  All  instruc¬ 
tions  are  the  same  as  for  card  number  four.  The  response 
functions  used  in  the  input  library  are  presented  in  Table 

IV.  The  functions  were  obtained  from  the  SMAUG  code. 

Card  No.  7: 

This  card  contains  one  number  in  14  format.  This  number 
is  the  number  of  aircraft  being  evaluated  and  must  be  between 
one  and  100.  This  card  must  be  supplied  by  the  user. 

EXAMPLE:  bbb4 

Card  No.  8: 

This  may  be  one  or  more  cards  and  supplies  the  aircraft 
position.  The  positions  are  defined  in  an  X,Y,Z  coordinate 
system  (in  that  order)  and  the  units  on  the  coordinates  are 
kilometers.  Each  card  contains  the  coordinates  for  two  air¬ 
craft  in  1X.1P6E11.4  format  and  will  normally  have  six 
numbers.  A  card  could  have  three  numbers  (one  aircraft 
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y^r-s^f-7 


Table  IV 


Gamma  Response  Functions 


Gamma  Group 

For  Silicon  Dose 
in  rads 

For  Tissue  Dose 
in  rads 

1 

2.8E-09  (a) 

2.2576E-09 

2 

2.28E-09 

1 . 937  IE-09 

3 

1.83E-09 

1.6436E-09 

4 

1.48E-09 

1 . 3901E-09 

5 

1.2E-09 

1 . 1812E-09 

6 

9.85E-10 

1 . 0102E-09 

7 

8.4E-10 

8.8378E-10 

8 

7.12E-10 

7 . 6726E- 10 

9 

6.1E-10 

6.6473E-10 

10 

5.05E-10 

5.5115E-10 

11 

4 . IE- 10 

4 . 4620E- 10 

12 

2.7E-10 

3.5710E-10 

13 

2.37E-10 

2 . 6009E-  10 

14 

1.65E-10 

1 . 7915E-  10 

15 

1.17E-10 

1 . 2257E- 10 

16 

7.25E-11 

6.2561E-11 

17 

9.75E-11 

3.2027E-11 

18 

4.15E-10 

4.7209E-11 

(a)  read  as  2.8  x  10 
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position)  if  only  one  aircraft  is  listed  on  card  number  seven 
or  if  the  card  is  the  last  card  for  an  odd  number  of  air¬ 
craft.  The  number  of  aircraft  positions  must  agree  with  the 
number  of  aircraft  specified  on  card  number  seven.  This  card 
(or  cards)  must  be  supplied  by  the  user. 

EXAMPLE: 

bbl.0000E+00bl.0000E+00bl. 0000E+0 lbO . OOOOE+OOb5 . OOOOE-Olb 1 .5000E  +  00 
bbl.2500E-01bi.OOOOE  +  OOb2.0000E+C)lbO.OOOOE+OObO.OOOOE+OOb2.3000E+00 
Card  No.  9: 

This  card  contains  one  or  two  numbers  in  214  format.  The 
number  specifies  the  units  of  aircraft  vulnerability  and  must 
be  a  number  from  one  to  eight.  If  the  first  number  is  five 
or  eight  no  second  number  is  required.  The  first  number  gives 
the  units  of  neutron  vulnerability  and  the  second  number  gives 
the  units  of  gamma  vulnerability.  The  meaning  of  the  numbers 
follow : 


Number 

1 

2 

3 

4 

5 

6 

7 

8 


Meaning 

2 

Total  neutron  fluence  (neutrons/cm  ) 

2 

Total  gamma  fluence  (gammas/cm  ) 
Neutron  tissue  dose  (rads) 

Gamma  tissue  dose  (rads) 
Neutron+gamma  tissue  dose  (rads) 
Neutron  silicon  dose  (rads; 

Gamma  silicon  dose  (rads) 
Neutron+gamma  silicon  dose  (rads) 


This  card  must,  be  supplied  by  the  user. 
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EXAMPLE:  bbblbbb2 

Card  No.  10: 

This  card  contains  one  or  two  numbers  in  1X,3P2E11.4 
format.  The  numbers  specify  the  aircraft  vulnerability. 

The  first  number  is  the  aircraft  neutron  vulnerability  and 
the  second  number  is  the  aircraft  gamma  vulnerability.  The 
units  must  be  those  specified  on  card  number  nine.  If  the 
first  number  on  card  nine  is  a  five  or  eight  only  one  number 
is  required  on  this  card.  This  card  must  be  supplied  by 
the  user. 

EXAMPLE:  bb 1 . 0000E  + lObl . 0000E+ 1 0 

Card  No.  11: 

.This  card  contains  four  numbers  in  1X,1P4E11.4  format. 
The  first  number  is  the  yield  of  the  nuclear  weapon  in  kilo- 
tons.  The  next  three  numbers  give  the  position  of  the  burst 
in  X  ,  Y ,  Z  coordinates  (in  that  order).  The  units  of  the 
coordinates  are  kilometers.  This  card  must  be  supplied  by 
the  user. 

EXAMPLE :  bbl . 0000E  +  0 lbl . 2500E-0 lbl . OOOOE-O lb  2 . 5C00E  +  0 1 
Card  No.  12: 

This  card  may  contain  two  numbers  in  214  format  and 
specifies  the  type  of  output  desired.  The  first  number  may 
be  blank  (or  zero),  one,  two,  or  three.  This  number  speci¬ 
fies  the  type  of  plot  output.  A  blank  means  no  plots  will 
be  output.  A  one  means  the  vulnerability  isofluence  or 
isodose  line  will  be  plotted  and  the  aircraft  will  be  located 
on  the  plot.  A  two  means  only  isofluence  or  isodose  lines 
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will  be  plotted.  A  three  means  both  the  one  and  two  options 
will  be  plotted.  The  second  number  may  be  blank  (zero)  or 
one.  This  number  specifies  the  type  of  printed  output.  A 
blank  is  the  normal  mode  and  will  result  in  a  short  printout 
giving  the  location  of  the  burst,  location  of  the  aircraft, 
and  if  the  aircraft  survived  with  gamma  and  neutron  levels 
experienced  by  the  aircraft.  A  one  will  give  the  same  out¬ 
put  as  the  blank  plus  the  detailed  group  by  group  fluence 
at  every  mesh  point.  The  one  option  generates  several  hundred 
pages  and  should  not  be  used  unless  the  group  fluences  are 
needed.  This  card  must  be  supplied  by  the  user. 

EXAMPLE:  bbb3bbbb 
Card. No.  13: 

This  card  may  contain  one  to  three  numbers  in  314  format. 
The  first  number  must  be  a  one  or  two  and  specifies  the  type 
of  nuclear  weapon.  A  one  is  a  fission  weapon  and  a  two  is  a 
thermonuclear  weapon.  The  second  number  may  be  a  blank  or 
a  one  and  specifies  the  source  of  the  weapon  output  spectrum 
for  neutrons.  If  this  number  is  blank,  an  unclassified 
default  spectrum  will  be  used.  This  default  spectrum  is 
contained  within  the  code  and  will  be  a  fission  or  thermo¬ 
nuclear  spectrum  depending  upon  the  first  number  on  this 
card.  Table  V  lists  the  default  neutron  spectra.  If  this 
number  is  one,  the  default  spectrum  will  not  be  used  and  a 
user  supplied  spectrum  will  be  used.  The  third  number  may 
be  a  blank  or  a  one  and  specifies  the  source  of  the  weapon 
output  spectrum  for  gammas.  If  this  number  is  blank,  an 
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Table  V 


Weapon  Neutron  Spectra 


Neutron  Group 


Fission  Spectrum 
(neutrons /ki loton) 

Thermonuclear  Spectrum 
(neutrons /ki lot on) 

3.92E+19  (a) 

6.001E+22 

2.233E+20 

2 . 176E+22 

8.7E+20 

1 . 1985E+22 

3.48E+21 

1 . 2495E+22 

8 . 705E+21 

1 . 4875E+22 

8 . 705E+21 

1 . 4875E+22 

1 . 49S1E+22 

1 . 4167E+22 

1 . 49S1E+22 

1 . 4167E+22 

1.4951E+22 

3 . 825E+22 

4.23E+22 

3 . 825E+22 

4.23E+22 

3 . 825E+22 

4 . 2325E+22 

7.9475E+22 

4. 2325E+22 

7 . 947SE  +  22 

3. 875E+21 

3. 1025E  +  23 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(a)  read  as  3.92  x  10 


GNE/PH/72-8 


y 


unclassified  default  spectrum  will  be  used.  This  default 
spectrum  is  contained  within  the  code  and  will  be  a  fission 
or  thermonuclear  spectrum,  depending  upon  the  first  number 
on  this  card.  Table  VI  lists  the  default  gamma  spectrum. 

If  this  number  is  one,  the  default  spectrum  will  not  be  used 
and  a  user  supplied  spectrum  will  be  used.  The  default 
spectra  are  from  the  SMAUG  computer  code.  If  both  the  second 
and  third  numbers  are  blank,  this  is  the  last  card  required. 
This  card  must  be  supplied  by  the  user. 

Card  No.  14: 

If  the  second  number  on  card  13  was  blank,  omit  this 
card.  If  that  number  was  a  one,  this  card  (actually  a  small 
deck)  will  contain  the  neutron  output  spectrum  of  the  weapon 
in  1X,1P7E11.4  format.  The  number  of  entries  must  agree 
with  the  number  of  neutron  groups  specified  on  card  one. 

The  group  structure  must  be  arranged  in  order  of  descending 
energy. 

Card  No.  15: 

If  the  third  number  on  card  13  was  blank,  omit  this 
card.  If  that  number  was  a  one,  this  card  (or  several  cards) 
will  contain  the  gamma  output  spectrum  of  the  weapon  in 
1X,1P7E11.4  format.  The  number  of  entries  must  agree  with 
the  number  of  gamma  groups  specified  on  card  number  one. 

The  group  structure  must  be  arranged  in  order  of  descending 
energy . 

Card  15  is  the  last  data  card  that  may  be  supplied  by 
the  user.  As  can  be  noted  in  the  above  description  of  the 
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Table  VI 


Weapon  Gamma  Spectra 


.  Group 

Fission  Spectrum 
( gammas /ki lot on) 

- 1 

Thermonuclear  Spectrum 
(gammas /kilo ton) 

1 

1 . 264SE+19  (a) 

6 . 3239E+18 

2 

5. 9019E+19 

2.S508E+19 

3 

1 . 0539E+20 

5 . 2693E+19 

4 

4 . 7555E  +  20 

2 . 3778E+20 

5  • 

4 . 7556E+20 

2 . 3778E+20 

6 

1 . 0718E+21 

5.3595E+20 

7 

1 . 0719E  +  21 

5 . 3595E+20 

8 

2 . 3562E+21 

1.I781E+21 

9 

3. 22E+21 

1 .615E+21 

10 

4 . 0838E  +  21 

2 .0419E+21 

11 

3. 7741E+21 

1 . 887E+21 

12 

4.512E+21 

2.2559E+21 

13 

S.2499E+21 

2.6249E+21 

14 

2. 2981E+21 

1 . 149E+21 

15 

2.2981E+21 

1.149E+21 

16 

2 . 7062E+21 

1 . 3531E+21 

17 

0 

0 

18 

0 

0 

(a)  read  as  1.2648  x  10*9 
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cards,  this  code  is  flexible  and  will  allow  the  input  library 
and  contained  data  to  be  replaced.  However,  if  any  of  the 
options  of  replacing  cross  sections,  response  functions,  or 
weapon  spectrum  are  exercised,  the  user  should  note  the 
following  caution.  These  three  selections  of  data  are  self- 
consistant  in  number  of  groups  and  the  energy  range  in  each 
group.  If  any  one  (or  more)  are  replaced,  this  self- 
consistency  must  be  retained.  That  is,  the  user  is  respon¬ 
sible  for  insuring  the  number  of  groups,  and  the  energy 
ranges  in  each  group,  agree. 

End  of  Job  Cards 

For  the  Wright-Patterson  Air  Force  Base  CDC  6600,  the 
user  must  supply  two  more  cards  following  the  data  card. 

The  first  is  a  multipunched  7/8/9  in  column  one.  The  second 
is  a  multipunched  6/7/8/9  4 n  column  one.  This  last  card  is 
orange.  Other  computers  will  require  different  end  of  job 
cards . 
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IV.  Results 

This  chapter  presents  the  results  obtained  from  the 
sample  problem  illustrated  in  Chapter  II/.  These  results 
are  compared  with  results  obtained  from  a  constant  density 
air  model  with  the  same  input  data.  The  results  are  also 
compared  to  those  given  in  the  Quick-Look  Radiation  charts 
(Ref  1).  These  charts  were  derived  from  SMAUG;  therefore, 
this  comparison  is  effectively  between  the  author's  code 
and  SMAUG. 

The  Problem 

The  problem  presented  in  Chapter  III  (as  sample  entries) 
is  reviewed  in  this  paragraph. 

Weapon  Parameters 

Position:  X  =  0.125  km,  Y  =  0.1  km,  Z  =  25  km 

Yield:  10  kilotons 

Type:  Fission 

Neutron  spectrum:  Input  library  (see  Table  V) 
Gamma  spectrum:  Input  library  (see  Table  VI) 
Aircraft  Parameters 
Number:  Four 

Positions:  No.  1,  X  =  1  km,  Y  =  1  km, 

Z  =  10  km 

No.  2,  X  =  0  km,  Y  =  0.5  km, 

Z  =  1.5  km 

No.  3,  X  =  0.125  km,  Y  =  1  km, 

Z  =  20  km 

No.  4,  X  =  0  km,  Y  =  0  km, 

Z  =  2.3  km 
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10  2 

Neutron  vulnerability:  10  neutrons/cm 

'  10  2 

Gamma  vulnerability:  10  gammas/cm 
Air  Parameters 

Cross  sections:  the  22  group  neutron,  18 

group  coupled  gamma  library 
(see  Appendix  C) 

Response  functions:  Input  library  (see 

Tables  III  and  IV) 


The  Output  from  the  Code 

The  output  required  was  the  short  printed  output  which 

states  what  happened  to  the  aircraft  and  four  plots;  two 

showing  neutron  and  gamma  isofluence  lines,  and  two  showing 

the  neutron  and  the  gamma  isovulnerability  line.  The 

printed  output  stated  that  aircraft  number  one  survived  and 

6  2 

experienced  a  neutron  fluence  of  4.8  x  10  neutrons/cm  and 

6  2 

a  gamma  fluence  of  5.4  x  10  gammas/cm  .  Aircraft  numbers 

two  and  four  survived  and  experienced  zero  neutron  and  gamma 

fluences.  Aircraft  number  three  was  killed  by  neutrons 

12  2 

with  a  fluence  of  2.2  x  10  neutrons/cm  .  The  four  plots 
are  presented  as  Figs.  2  through  5. 


Comparison  with  Resul  ts  from  a^  Constant  Density  Air  Model 

Another  air  model  was  selected  for  comparison  purposes. 
This  model  was  a  constant  density,  homogeneous  composition, 
infinite  atmosphere  model.  The  same  input  data  was  used. 

A  program  was  written  for  this  model  with  the  same  require¬ 
ment  on  output.  The  printed  output  stated  all  four  aircraft 
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FIG.  2  NEUTRON  ISC FLUENCE  LINES  ( N/CM2  ) 


EXPONENTIAL  AIR 
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were  killed.  Aircraft  number  one  experienced  a  neutron 

11  2 

fluence  of  6.6  x  10  neutrons/cm  ,  aircraft  number  two 

11  2 

experienced  a  neutron  fluence  of  1.9  x  10  neutrons/cm  , 

aircraft  number  three  experienced  a  neutron  fluence  of 
12  2 

3.1  x  10  neutrons/cm  ,  and  aircraft  number  four  experi- 

11  2 

cnced  a  neutron  fluence  of  2.1  x  10  neutrons/cm  .  The 
four  plots  are  presented  as  Figs.  6  through  9. 

The  neutron  fluences  calculated  at  the  four  aircraft 
positions  are  greater  for  the  constant  density  model.  This 
difference,  however,  is  expected  since,  in  the  actual 
atmosphere,  air  density  is  decreasing  with  altitude.  There¬ 
fore,  in  the  actual  atmosphere,  any  fixed  position  below  the 
burst  should  experience  a  lower  fluence  than  that  predicted 
for  a  constant  density  model.  The  exponential  air  model 
selected  for  the  code  is  a  reasonable  approximation  to  the 
actual  atmosphere,  therefore,  the  fluences  calculated  with 
this  code  for  fixed  positions  below  the  burst  should  be 
lower  than  those  predicted  for  the  constant  air  model. 

Two  observations  can  be  made  by  examining  the  plots 
from  the  two  codes.  Figure  6,  the  neutron  isofluence  lines 
for  the  constant  density  model,  and  Fig.  2,  the  neutron 
isofluence  lines  for  the  exponential  air  model,  show  the 
effect  of  the  atmospheric  model  on  neutron  fluence.  The 
isofluence  lines  for  the  constant  density  model  are  circles. 
The  isofluence  lines  for  the  exponential  air  model  are  not 
circles.  These  lines  are  close  together  in  the  region  of 
greater  air  density  and  diverge  as  the  air  density  decreases. 
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LINE  C  N/CM2 )  CONSTANT  0EN5ITY  AIR 
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This  observation  can  also  be  made  by  coirparing  Figs.  4  and 

8,  the  gamma  isofluence  lines. 

The  second  observation  that  can  be  made  is  that  the 

isofluence  lines  at  a  height  coaltitude  with  the  burst  differ 

between  the  two  models.  This  difference  can  be  easily  seen 

by  comparing  Figs.  3  and  7,  the  neutron  fluence  vulnerability 

line,  or  by  comparing  Figs.  5  and  9,  the  gamma  fluence 

vulnerability  line.  From  Fig.  3,  the  vulnerability  level 
10  2 

of  10  neutrons/cm  is  at  about  18  km  coaltitude  for 

exponential  air.  From  Fig.  7,  the  vulnerability  level  of 
10  2 

10  neutrons/cm  is  at  about  44  km.  Figures  5  and  9  show 

similar  results  for  the  gamma  vulnerability  line  of 
10  2 

10  .gammas/cm  .  The  coaltitude  distance  for  exponential 
air  is  about  12  km  but  for  constant  density  air  is  about 
25  km.  Therefore,  in  this  example,  the  constant  density 
air  model  overpredicts  the  coaltitude  distance  by  a  factor 
of  two. 

Comparison  with  the  Quick-Look  Radiation  Charts 

Neutron  and  gamma  fluences  at  the  aircraft  positions 
were  also  calculated  using  the  Quick-Look  Radiation  Charts 
(Ref  1).  These  charts  were  derived  from  SMAUG.  The 
expected  values  from  these  charts  should  be  between  the 
exponential  air  model  and  the  constant  density  air  model 
because  the  SMAUG  data  base  is  from  a  constant  density  model. 
SMAUG  compensates  for  the  actual  change  in  air  density  by 
mass  integral  scaling  along  the  line  of  sight  between  the 
burst  and  the  receiver. 


65 


GNE/PH/72-8 


g 

The  fluence  values  from  these  charts  were  1.6  x  10 
2  8  2 

neutrons/cm  and  3.8  x  10  gammas/cm  for  aircraft  number 

12  2  11  2 
one,  3  x  10  neutrons/cm  and  3  x  10  gammas/cm  for 

aircraft  number  three.  Values  for  aircraft  numbers  two  and 

four  could  not  be  obtained  because  they  were  out  of  the 

range  of  the  charts.  The  results  of  all  these  calculations 

are  summarized  in  Tables  VII  and  VIII. 


Table  VII 

Neutron  Fluences  at  the  Aircraft  Positions 
as  Calculated  by  the  Three  Models 


Aircraft 

Number 

2 

Neutron  Fluence  (neutrons/cm  ) 

Exponential  Air 

Quick- Look 
Charts 

Constant  Density 
Air 

1 

vO 

O 

rH 

X 

CO 

• 

rf 

1.6  x  108 

6.6  x  10U 

2 

0 

(a) 

1.9  x  1011 

3 

12 

2.2  x  101Z 

3  x  1012 

12 

3.1  x  10A 

4 

0 

(a) 

2.1  x  1011 

(a)  Out  of  the  range  of  the  charts 


The  comparison  between  the  three  calculations  do  show 
the  Quick-Look  results  between  those  of  exponential  air  and 
constant  density  air. 

Validity  of  Results 

The  general  appearance  of  the  results  seem  valid.  The 
results  show  the  expected  effect  of  the  decreasing  air 
density  with  height.  Furthermore,  when  these  results  are 
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Table  VIII 


Gamma  Fluences  at  the  Aircraft  Positions 
as  Calculated  by  the  Three  Models 


Aircraft 

Number 

2 

Gamma  Fluence  (gammas/cm  ) 

Exponential  AiT 

Quick-Look 

Charts 

Constant  Density 
Air 

1 

5.4  x  106 

3.8  x  108 

8  x  1010  (c) 

2 

0 

(b) 

3  x  1010  (c) 

3 

1011  (a) 

3  x  1011 

>  1011  (c) 

4 

0 

(b) 

3  x  1010  (c) 

(a)  Not  calculated,  estimated  from  Fig.  4 

(b)  Out  of  range  of  the  charts 

(c)  •  Not  calculated,  estimated  from  Fig.  8 


compared  to  the  constant  density  air  model  and  to  the 
Quick-Look  charts,  they  appear  to  be  correctly  predicting 
increased  (or  decreased)  fluences  where  applicable.  The 
actual  validity  of  the  numbers  can  not  be  determined  since 
no  experimental  observations  exist  for  these  type  of  calcula¬ 
tions.  Some  checks,  other  than  comparison  with  experimental 
data,  can  be  made  to  estimate  the  validity  of  the  results. 
First,  a  check  for  conservation  of  particles  can  be  made. 

The  conservation  relationship  is  that  the  number  of  particles 
entering  a  specified  volume  plus  the  number  of  particles 
gained  inside  the  volume  is  equal  to  the  number  of  particles 
leaving  the  volume  plus  the  number  of  particles  lost  inside 
the  volume .  This  check,  once  made,  would  establish  that 
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the  difference  equations  used  in  this  report  do  conserve 
particles  and  are  therefore  valid.  Another  check  that  can 
be  made  is  to  let  the  atmospheric  model  used  approach  a 
constant  density  model.  This  can  be  done. by  changing  the 
scale  height  of  the  atmosphere  to  a  very  large  number.  If 
the  same  example  used  in  this  report  is  then  run,  the 
results  should  approach  those  obtained  from  a  constant 
density  model. 
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V.  Conclusions  and  Recommendations 


Conclusions 

Definite  conclusions  can  not  be  made  since  only  one 
sample  problem  has  been  run  on  the  computer  code.  One 
tentative  conclusion,  however,  can  be  made  subject  to  con¬ 
firmation  by  more  computer  runs.  The  codes  that  consider 
coaltitude  burst  and  receiver,  if  based  in  the  constant 
density  air  assumption,  overestimate  the  gamma  and  neutron 
fluence.  The  reason  for  the  overestimation  is  most  likely 
the  implicit  assumption  that  scatter  from  above  is  equal  to 
scatter  from  below  in  the  constant  density  air  model.  In 
actuality,  due  to  the  density  variation,  scatter  from  above 
would  be  less  than  scatter  from  below.  Therefore,  the 
exponential  model  correctly  indicates  lower  coaltitude 
f luences . 

Another  conclusion  that  can  be  drawn  is  that  the 


const 

ant  density  air 
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odel  poorly 

estim 

ate 
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gamma 
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fully  satisfied.  The  running  time  of  the  code,  20  minutes 
of  central  processor  time  plus  one  hour  of  input/output 
time,  does  not  compete  favorably  with  Sf'AUG’s  running  time 
of  several  seconds,  but  it  is  much  faster  than  the  Monte 
Carlo  codes  mentioned  in  Straker's  report  (Ref  5). 

Recommendations 

Based  on  the  above  conclusions,  the  author  has  five 
recommendations.  The  first  is  that  a  conservation  check  be 
performed,  using  the  geometry  developed,  to  insure  the 
validity  of  the  difference  equations.  The  second  is  that 
the  code  be  rerun  with  an  increased  scale  height  to  determine 
if  the  results  do  approach  those  obtained  with  a  constant 
density  air  model.  The  third  recommendation  is  that,  if  the 
above  checks  are  successfully  made,  the  code  be  rerun  at 
heights  of  burst  for  which  Monte  Carlo  data  exists  and  a 
comparison  made  between  the  Monte  Carlo  results  and  those 
obtained  from  this  code.  The  fourth  recommendation  is 
that,  if  the  above  three  checks  still  show  the  code  to  be 
valid,  the  code  be  given  to  an  experienced  programmer  to 
revise  in  order  to  decrease  the  run  time.  The  final 
recommendation  is  to  then  run  the  revised  code  at  a  series 
of  burst  heights  in  order  to  generate  a  broad  data  base. 

This  data  base  can  then  be  used  to  write  a  SMAUG-like 
program  which  should  fully  satisfy  the  requirement  for  an 
inexpensive  (quick  running)  code. 
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Appendix  A 

An  Aigorithn  for  the  Solution  of  Matrix  Equations 
with  Block  Tri-Piagonal  Matrix  of  Coefficients 

A  natrix  equation  of  the  forn 

AF  =  S  (1) 

where  A  is  an  M  x  M  natrix,  F  is  an  M  x  1  natrix,  and  S  is 
an  M  x  1  natrix,  can  be  directly  solved  for  £  by  finding  A 
inverse  and  nultiplying  £  by  A  inverse.  The  solution  is 
therefore 

£  =  A_1S  (2) 

This  is  the  nornal  nethcd  of  solution  used  in  nost  conputer 
prograss.  However,  when  N  is  large,  this  cethod  is  imprac¬ 
tical.  Two  reasons  for  the  inpracticality  of  this  nethod 
are  the  large  conputer  core  and  the  long  conputer  run  tine 
required.  For  ezanple,  if  N  is  100 0,  a  computer  CGre  of 
1,000,000  words  will  be  required  just  to  store  the  inverse. 
The  core  requirement  can  be  reduced  by  using  magnetic  tapes 
fr *■  storage;  however,  the  use  of  magnetic  tapes  increases 
the  conputer  run  tine.  Xunerous  algorithms  have  therefore 
been  developed  to  avoid  the  problems  of  solution  by  direct 
inversion. 

These  algorithms  depend  os  the  composition  of  the 
coefficient  matrix  A.  If  A  is  a  block  tri-diagonal  matrix, 
an  algorithm  has  been  described  by  Winchester  (Ref  8:57-61) 
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that  greatly  reduces  the  coaputer  core  size  and  run  tiae  to 
solve  Eq  (1).  This  algoritha  is  described  in  the  following 
paragraphs . 

If  A  is  block  tri-diagonal,  it  can  be  partitioned  such 


that 


where  DIA^,  UP.  and  LOST.  are  all  square  aatrices. 
also  be  partitioned  such  that 


F  and  S  can 


(4) 


S  = 


(5) 


GNE/PH/72-8 


where  SM.  and  FM.  are  also  natriccs 
x  1 

Next,  factor  A  into  natrices  KM  and  QM  such  that 

KM  QM  =  A  (6) 

where 


Carrying  out  the  oatrix  nultiplication  of  Eq  (6)  and 
equating  conponents  to  corresponding  diagonal  conponents  of 
A  yields  the  following  relations: 

Kj  =  DIAJ  (9) 


K.  =  DIA.  -  LOb.  q.  ,  (10) 

ill  l-l 
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By  equating  components  with  corresponding  off-diagonal 
components  of  A  the  following  relation  is  also  obtained: 


q.  =  (K-)'1UP. 

ni  i  x 


(11) 


Next,  define  *»  matrix  GM 


such  that 


GM  = 


/m 


\GN 


(12) 


WM  GM  =  S 


(13) 


Therefore 


GM  =  QH  F 


(14) 


Carrying  out  the  matrix  multiplication  of  F.q  (13) 
and  equating  the  components  with  corresponding  components 
of  S  yields  the  following  relations: 


Gi  =  <V  1  SM1 


(IS) 


G.  =  (K.)_1(SM.  -  LOtf.G.  t)  1  <  i  <  N  (16) 
l  l  '  l  x  l-l 


The  same  operation  with  Eq  (14)  yields 


F\ 


FH. 

f 


S 


G. 

l 


-  q.FM.  . 
U  l*: 


1  <  i  <  li 


(17) 

(IS) 
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Equations  (14)  and  (15)  are  the  equations  defining  the 
unknown  na.tr ix  F. 

The  entire  algorithm  to  solve  for  F  is  listed  below. 

(1)  Start  with  i  =  1 

(2)  Solve  for  W. 


W1  =  DIAj 


W.  =  DIA.  -  LON. q.  . 

l  l  iMi-1 


1  <  i  <  N 


(3)  Solve  for  q. 

°*i  =  (Wi^*luPi 

(4)  Solve  for  G.^ 

G1  =  C»1)‘1SM1 

6.  =  (W.)"1(SM.  -  LOW . G .  ,)  1  <  i  <  N 

(5)  Return  to  Step  (2)  with  next  i  and  repeat  until 
ail  conponents  are  calculated. 

(6)  After  all  the  q.  matrices  and  the  G.  natrices  have 
been  calculated,  solve  for  the  FM.  components  of  the  unknown 
matrix  F.  Start  with  i  =  N. 


F!,S  -  S 

F«i  =  G.  -  ,  FH.M 


1  <  i  <  N 
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Appendix  B 


Description  and  Listing  of  the  Code 


Introduction  .  .  . 
Program  MAIN  .  .  . 
Subroutine  CROSS  . 
Subroutine  DOSE 
Subroutine  AIRCR  . 
Subroutine  GAMNEUT 
Subroutine  SELECT 
Subroutine  MESH 
Subroutine  LOCATE 
Subroutine  VIRGIN 
Subroutine  ADD  .  . 
Subroutine  OUT  .  . 
Subroutine  CONVERT 
Subroutine  CHECK  . 
Subroutine  NOTICE 
Subroutine  MAP  .  . 
Subroutine  CONTOUR 
Subroutine  SETUP  . 


Page 

78 

78 

81 

83 

85 

87 

98 

101 

104 

107 

110 

112 

114 

117 

120 

123 

128 

131 


BLOCK  DATA 
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Introduction 

This  appendix  describes  and  presents  the  code.  Each 
subroutine  is  described  separately  with  a  glossary  and  listing. 
Some  program  modifications  that  can  be  made  with  the  exchange 
or  addition  of  a  few  cards  are  also  discussed.  The  code  is 
written  for  the  Wright-Patterson  Air  Force  Base  CDC  6600  in 
FORTRAN  Extended.  In  addition,  the  code,  as  written,  should 
only  be  run  on  a  1700  terminal  of  the  CDC  6600  that  has 
on-line  plot  capability.  The  program  can  easily  be  converted 
to  run  at  the  computer  center,  rather  than  at  a  terminal, 
with  a  few  easy  modifications.  These  modifications  are  dis¬ 
cussed  in  subroutine  MAP.  If  the  code  is  run  on  another 
CDC  computer,  the  library  matrix  and  plot  subroutines  used 
by  the  code  may  be  different  and  require  modifications  to 
the  code. 

Program  MAIN 

This  is  the  main  program  and  is  a  substitute  for  the 
master  nuclear  effects  program  developed  by  Capt  DeRaad 
(Ref  2).  The  program  is  responsible  for  feeding  the  cross 
section  data,  response  function  data,  aircraft  data,  weapon 
yield  and  position,  and  output  options  to  subroutine  GAMNEUT. 
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FORTRAN 

ACPOS 

ARRAY 

BURST 

LMAP 

LOUT 

MODE 

NGAM 

NGROUP 

NNEUT 

NUMBER 

STMAX 

STMIN 

VUL 

YIELD 


Glossary: 

An  array  that  stores  the  aircraft  positions 
in  X,  Y,  Z  coordinates.  Units  are  kilometers 

An  array  that  specifies  a  live  aircraft  or 
one  killed.  A  one  means  the  aircraft  is  live 
while  a  zero  means  the  aircraft  has  been 
killed. 

An  array  that  store  the  burst  position  in  X, 
Y,  Z  coordinates.  Units  are  kilometers. 

A  numbeT  indicating  the  p;.ot  option. 

A  number  indicating  the  printed  output  option 

An  array  indicating  the  units  of  the  figures 
given  for  aircraft  vulnerability. 

The  number  of  gamma  groups. 

The  number  of  gamma  plus  neutron  groups. 

The  number  of  neutron  groups. 

The  starting  number  of  aircraft. 

The  maximum  macroscopic  total  air  cross 
section. 

The  minimum  macroscopic  total  air  cross 
section. 

An  array  that  stores  the  aircraft  neutron  and 
gamma  vulnerability. 

The  weapon  yield  in  kilotons. 
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PROGRAM  MAIN(INFUT,TAFE9,TAPEi: ,TAPEil, OUTPUT, PLOT) 
INTEGER  ARRAY ( i  j  G ) 

DIMENSION  AC°03 (3 ,ir 3 > , BURST ( 3) , VUL (2 > , MODE (2  ) 
COMMON  I 

READ  IN  T»£  NUMPER  OF  GAMMA  AND  NEUTRON 
GROUPS.  READ  IN  THE  GROUP  CROSS  SECTIONS. 

READ  IN  NEUTRON  AND  GA*MA  RESFCfSE  FUNCTTCNS. 

CALL  CROSS  <STMAX,STMIM,NNEUT,NGAM,NGRCUP) 

CALL  DOS.E(NNEUT,NGAM) 

RFAD  IN  AIRCRAFT  DATA,  THE  NUMBER  CF  AIRCRAFT, 
POSITIONS,  VULNFR A °ILITY .  ALSO  SET  UP  AN 
APRAY  TO  KE-e  TRACK  OF  AIRCPA-T  1  HAT  ARE 
KILLED. 

CALL  AIRC°( ARRAY, ACPOS, VUL, MODE, NLMEER) 

READ  IN  DURST  LOCATION  AND  YTELC. 

READ  130, YIELD, (8UPST(I) ,1=1,3) 

PEAD  IN  PLOT  ANO  OUTPUT  OPTIONS. 

LMAP=C  (3LANK)  MEANS  NO  PLOT. 

LMAPsi  MEANS  A  V'JLNEFASILITY  ISCFLUENCE  OR 
DOSE  LINE  HILL  *E  PLOTTED  ANC  THE 
AIRCRAFT  LOCATED  ON  THE  PLOT. 

LMAP=2  MEANS  THE  ISOFLUENCE  LINES  HILL 
8E  PLOTTED. 

LMAP=3  MSANS  DOTH  LMAF=1  AND  LMAF=2  OPTIONS 
HILL  PE  PLOTTED. 

LOUT="  (BLANK)  MEANS  SHORT  (NORMAL)  OUTPUT. 
LOUT=l  MEANS  A  DETAILED  CUTPUT. 

RRAD  l.,i,LMAP,LCUT 


CALCULATE  THE  NEUTRCN  ANC  GAMMA  E NVIRCMEN'T. 

CALL  GAMNEUT  (MjMorp, ACFOS, ARRAY, DURST, VUL, MOD E,NNEUT, NGAM, STMAX, ST 
AMIN, LMAP, LOUT, N3RCUF, YIELC) 

1^0  FC°MfiT(iX,iPoEU.<.) 

131  FOPMAT (1314) 

STOP 

END 
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Subroutine  CROSS 

This  subroutine  reads  in  the  sea  level  macroscopic 
cross  sections,  determines  the  minimum  and  maximum  total 
cross  sections,  and  stores  all  cross  sections  on  a  disk 
file.  This  subroutine  should  be  included  in  the  master 
nuclear  effects  program. 


FORTRAN  Glossary: 

NGAM  The  number  of  gamma  groups. 

NGROUP  The  number  of  gamma  plus  neutron  groups. 

NNEUT  The  number  of  neutron  groups. 

STMAX  The  maximum  macroscopic  total  air  cross 

section . 


STM  IN 


The  minimum  macroscopic  total  air  cross 
section. 


XSECT 


An  array  containing  the  neutron  and  gamma 
sea  level  air  macroscopic  cross  sections. 
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SUBROUTINE  CROSS (STMAX ,STMIN ,NNEUT, NGAM , NGROUc) 

\ 

THIS  SUBROUTINE  REACS  IN  THE  GAMMA  AN? 
NEUTRON  C*CSS  SECTIONS  AND  STORES'  THEM  ON 
TAPE  9. 


NGROUD  IS  THE  NUM°EP  CP  GAMMA  ANC  NEUTRON 
GROU°S • 


STMAX  IS  THE  LARGEST  TOTAL  CRCSS  SECTION. 

STMIN  IS  THE  SMALLEST  TOTALCROSS  SECTION. 

K  IS  THE  NUMBER  OF  CRCSS  SECTIONS  IN  ANY 
ONE  GROUP. 


IN  ANY  GROUP  G,  THE  CROSS  SPCTICNS  ARE 
.  ARRANGED  IN  THE  FOLLOWING  ORDER. 

1  IS  TRANSPORT 

2  IS  REMOVAL 

3  IS  TOTAL 

4  IS  SCATTER  (G  TC  G) 

5  IS  SCATTER  (G-i  TO  G> 

6  IS  SCATTER  (G-2  TO  G) 

7  IS  SCATTER  (G-3  TO  G) 

ETC. 


NOTE..  IF  OTHER  THAN  THE  SUPPLIED  DATA  IS 
USED,  THE  NUMBER  "F  GROUPS  AND  ENERGY  0 AN CS 
MUST  A G P E F.  WITH  TWAT  SUPPLTED  AS  INPUT  FOR 
DOSE  CALCULATIONS  AND  WEAPCN  SOLRCE  SPECTRUM. 

COMMON  XSECT(43,4') ,1, J,M 
REWIND  9 

READ  k.NNEUT.NOAM 
NGROUP=NNEUT+NGAM 
M=NGROUP+? 

READ  1,  ((XSECT(I,J),I=1,M),J=1,NGRCUF) 

DO  2  J?1,NGR3U° 

2  WRITE  (9)  (XSECT(I,J) ,I=1,M) 

STMAX=XSFCT (3,1) 

STMIM=XSECT (7,1) 

DO  3  T=l,NGROUP 

STMAX=AM AX1 (STMAX,XSE0T(3,I) ) 

?  STMTN=AMIN1  ( ST.-IIN  , XcEf'  T  C7  ,  I )  ) 

1  FORMAT (IX, 1P7E11.4) 

4  FORMAT(2IM 

END 
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Subroutine  DOSE 

This  subroutine  reads  in  the  response  functions  that 
are  used  to  convert  multigroup  fluence  to  silicon  dose  in 
rads  or  tissue  dose  in  rads.  These  response  functions  are 
then  stored  on  a  disk  file.  This  subroutine  should  be 
included  in  the  master  nuclear  effects  program. 

FORTRAN  Glossary: 

NGAM  The  number  of  gamma  groups. 

NNEUT  The  number  of  neutron  groups. 

SILGAM  An  array  containing  the  response  functions 

to  convert  group  gamma  fluences  to  silicon 
dose. 

SILNEUT  An  array  containing  the  response  functions 

to  convert  group  neutron  fluences  to  silicon 
dose. 

TISGAM  An  array  containing  the  response  functions 

to  convert  group  neutron  fluences  to  tissue 
dose. 

TISNEUT  An  array  containing  the  response  functions  to 

convert  group  neutron  fluences  to  tissue  dose. 
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SUBROUTINE  OOSE  ( NNEUT  ,  Nr, AH) 

THIS  SUBROUTINE  READS  IN  THE  NEUTRON  AND 
OAM 'I A  RESPONS-  FUNCTIONS  FOR  CONVERTING 
FLUENCE  TO  TISSUE  OOSE  OR  SILICON  COSE  IN 
PAOS. 

NOTE..  IF  OTwrp  THAN  THE  SUPFLIEC  CATA  IS 
USES,  THE  NUMBER  OF  GROUPS  AND  ENERGY  BANCS 
MUST  AGREE  WITH  THAT  SU°RLIEC  AS  INPUT  FOR 
GROUP  CROSS  STCTICNS  AND  WEAPON  SOURCE 
SPECTRUM. 

THE  RESPONSE  FUNCTIONS  ARR  STCREC  ON  TAPE  9. 

SILNEUT  IS  THE  RESPONSE  FUNCTION  FCR  CONVERTING 
NEUTRON  FLUENCF  TO  SILICON  COSE  IN  RADS. 

•  SI LG  AM  IS  THE  RESDONSE  FUNCTION  FCR  CONVERTING 
GAMMA  FLUENCE  TO  RILICON  ^OSE  IN  RADS. 

TISNE’JT  IS  THE  P£S°ONSE  FUNCTION  FOR  CCNVCPTING 
NEUTRON  FLUENT  TO  TISSUE  POSE  IN  RADS. 

TISGAM  IS  the  P £ S P C N S r  FUNCTICN  FCR  CONVERTING 
GAMMA  FLUENCE  TO  TISSUE  COSE  IN  RACS. 

COMMON  SILNEUT  (  22 )  ,SILGAM(  13)  ,TISNEUT  (22)  , TISGAM (13)  ,  I 
READ  1, (SILNEUT (I) ,T=1, NNEUT) 

WRITE  (9)  (SILNEUT(I) ,T  =  l, NNEUT) 

READ  1,  (SIL3AM(I),Isl,NGAM) 

WRITE  (9)  (SILGAM  (I) ,I  =  1,NGAM) 

RFAD  1,  (TISNEUT(I) ,I=1,NNCUT) 

WRITE  (9)  (TTSNEUT(I)  ,T  =  1,NNE!JT) 

READ  1,  (TISGAM(I), 1=1, NOAM) 

WRITE  (9)  (TISGAM(I),T=1,NGAM) 

REWIND  9 
1  F0RMAT(1Y, 1P7E11. A) 

RETURN 

END 
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Subroutine  AIRCR 

This  subroutine  reads  the  aircraft  data:  positions, 
vulnerability,  and  vulnerability  units.  It  also  sets  up 
the  array  that  keeps  track  of  whether  the  aircraft  survives 
or  is  killed.  This  subroutine  is  a  substitute  for  the  data 
that  would  be  obtained  from  the  master  nuclear  effects 
program. 

FORTRAN  Glossary: 

ACPOS  An  array  that  stores  the  aircraft  positions 

in  X,  Y,  Z  coordinates.  Units  are  kilometers. 

An  array  that  specifies  a  live  aircraft  or 
one  killed. 

An  array  indicating  the  units  of  the  figures 
given  for  aircraft  vulnerability. 

The  starting  number  of  aircraft. 

An  array  that  stores  the  aircraft  neutron  and 
gamma  vulnerability. 


ARRAY 

MODE 

NUMBER 

VUL 
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SUBROUTINE  AI°CR  (ARRAY, ACFCS, VUL ,MOCE , NUMBER) 
INTEGER  ARRAY (1 j3) 

DIMENSION  ACPOS{?,ir3) ,VUL (2) ,MODE (2) 

READ  1-11,  'JUMPER 
DO  i  I=i f NUMRER 
1  ARRAY (I) =1 

READ  iu }, ( (ACPOS (I, J) , T  =  i,3) ,J=1,NUM?ER> 

READ  ICi,  (M0'“!E(I)  ,1=1,2) 

READ  ltJf :VUL(I) ,1=1,2) 

1 0 Ci  FORMAT  <1X,1P6£11»A) 

101  FORMAT (101  A) 

END 
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Subroutine  GAMNEUT  ♦. 

This  subroutine  is  the  heart  of  the  code  and  is  the 
subroutine  designed  for  incorporation  into  the  master  nuclear 
effects  program.  This  subroutine  calculates  the  scattered 
neutron  and  gamna  fluences.  All  the  following  subroutines 
are  called  by  this  subroutine  and  have  functions  of  calcu¬ 
lating  the  preliminary  information  needed  for  the  fluence 
calculations  or  converting  the  data  to  a  form  suitable  for 


data  output 


Subroutine 

A 

ABAR 

ACPOS 

B 


DIA 


GAMNEUT  Glossary: 

A  coefficient  of  the . difference  equation 

A  coefficient  of  the  difference  equation 

An  array  that  stores  the  aircraft  positions 
in  X,  Y,  Z  coordinates 

A  coefficient  of  the  difference  equation 


& 

BBAR 

A  coefficient  of  the  difference  equation 

1 

BURST 

An  array  that  stores  the  burst  position  in 

h 

£ 

X,  Y,  Z  coordinates 

i 

tx 

fU 

C 

A  coefficient  of  the  difference  equation 

* 

s 

D 

The  diffusion  coefficient 

f. 

1 

DEIR 

The  radial  mesh  interval 

1 

1 

j* 

DELRIOW 

The  value  of  DELR  equivalent  to  the  value  of 

& 

1  3 

the  x  interval  in  the  highest  row  in  x  using 

% 

1 

the  nonorthogonal  coordinate  system 

I 

DELX1 

The  x*  interval 

fe¬ 

ll 

DELX3 

3 

The  x  interval 

A  packed  matrix.  This  is  one  of  the  sub¬ 
matrices  used  in  the  block  tri-diagonal 
algorithm 
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FLU 

G 

GROUP 

G1 

G2 

G3 

G4 

H 

HOB 

KCOL 

LMAP 

LOUT 

LOWER 

MODE 

NGAM 

NGROUP 

NHOR 

NNEUT 

NUMBER 

NUP 

POS 

Q 

S 


An  array  that  stores  the  calculated  group 
scattered  fluences 

One  of  the  matrices  used  in  the  block  tri¬ 
diagonal  algorithm 

A  counter.  Its  value  is  that  of  the  group 
for  which  scattered  fluence  is  being  calcu¬ 
lated 


Coefficients  of  the  difference  equation 
The  scale  height  of  the  atmosphere 
The  height  of  burst 

An  array  used  to  unpack  the  LOWER  matrix  for 
the  row  connecting  the  two  coordinate  systems 

A  number  indicating  the  plot  option 

A  number  indicating  the  printed  output  option 

A  packed  matrix.  This  is  one  of  the  sub¬ 
matrices  used  in  the  block  tri-diagonal 
algorithm 

An  array  indicating  the  units  of  the  figures 
given  for  aircraft  vulnerability 

The  number  of  gamma  groups 

The  number  of  gamma  plus  neutron  groups 

The  number  of  horizontal  mesh  points 

The  number  of  neutron  groups 

The  starting  number  of  aircraft 

The  number  of  vertical  mesh  points 

An  array  that  stores  the  aircraft  positions 
in  r,z  coordinates 

One  of  the  submatrices  used  in  the  block 
tri-diagonal  algorithm  . 

One  of  the  submatrices  used  in  the  block 
tri-diagonal  algorithm 
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SIGR 

SIGS 

SIGT 

SIGTR 

SOURCE 

SPECT 

CPPER 

XI 

X1MIN 

X3 

X3HOB 

X3MIN 

X3SW 

YIELD 


The  group  air  macroscopic  removal  cross 
section 

The  group  air  macroscopic  scatter  cross 
section 

The  group  air  macroscopic  total  cross  section 

The  group  air  macroscopic  transport  cross 
section 

The  total  number  of  neutrons  or  gammas  output 
from  the  weapon  for  a  group 

The  number  of  neutrons  or  gammas  per  kiloton 
output  from  the  weapon  for  a  group 

A  packed  matrix.  One  of  the  submatrices  used 
in  the  block  tri-diagonal  algorithm 

The  value  of  the  x*  coordinate 

The  minimum  value  of  the  x*  coordinate 

3 

The  value  of  the  x  coordinate 

The  height  of  burst  expressed  in  terms  of  x3 

The  minimum  value  of  x3 
3 

The  value  of  x  when  the  coordinate  systems 
are  switched 

The  weapon  yield  in  kilotons 
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subroutine  c amneut (numbep  accos, array, purst, vul,mooe,nne'j 

ASTMAX,  STMTN, LMAR,  LOUT, NGRi'UP, YIELD) 

this  SUBROUTINE  DETERMINES  THE  HULTIGRQUP 
GAMMA  AND  NEUTRCN  FLUENCE 

RRAL  LOWER (6 i , 3) 

INTEGER  ARRAY (13 J), GROUP 

OIMENSION  VUL(2)  ,  MODE  (?)  , ACPOS  (3,1j C )  ,EURST<  3),  SPECKS.  )  , 
AG)  ,P0R(2,i  JO 

COMMON  G(*jj,123)  ,C  (50,60  ,01 A  (6  3,3)  , UPPER (60, 3)  ,KC0L(6O  , 
AS<6-,120 
EPUIVALENCE(S,FLU) 

OFTERMINE  THE  WEAPON  OUTPUT  SPECTRUM. 

CALL  SELECT (SPECT,NNEUT,MGAM,NGROUP) 

LOCATE  HOB  (IN  UNITS  CF  CM.)  ANC  SET  U» 

SCALE  HEIGHT  H  AND  PI.  IF  HC°  IS  LESS  THAN 
7  KM.  OR  GREATER  THAN  100  KM.,  NO  CALCULATIONS 
WILL  BE  MADE. 

HOB=BURST (  2)  *1.0E+r5 

IF(H0B.LT.7.AE  +  35.0P.H0e.GE.l: J.0E  +  05)  301,302 
301  PRINT  5j1 
PRINT  501 

5 J Q  FORMAT (////, 4X,*wqp  is  LESS  THAN  7  KM.  OR  GREATER  THAN  1. 

501  FORMAT  (LX,  *NG  GAMMA  0°  NEUTRON  CALCULATIONS  W^ll  nE  MACE* 
GO  TO  Zu2 

332  H  =  7 . C  239E+  05 
RI=3.1415B27 
REWIND  11 

OFTERMINE  THE  MESH 

call  MESH(H,STMIN,STMAX,HCB,X1MIN,0ELX1,NH0R, L,X3MIN,CELX 
AX3H0R ) 

LOCATE  THE  AIRCRAFT  IN  THE  MESH.  IF  NO 
AIRCRAFT  ARE  IN  THE  MESH  A  Rr  A ,  ALL  AIRCRAFT 
ARE  ASSUMED  TO  HAVE  SLRVIVED.  NC  FURTHER 
CALCULATIONS  WILL  nE  MADE. 

CALL  LOCA TE (NUMBER, ACcOS, QURST,X 1MIN, CELX1,NHCP, H,X?M IN, I 
AL  ,POS  ,N.RET  .ARRAY) 

IF(M»ET.  £0.0303, 304 

333  PRINT  5n 2 
PRINT  5.1 

502  FORMAT<////,4X,*ALL  AIRCRAFT  APE  CUTSICE  THE  AREA  OF  C-AM* 
AUTRONS*) 

GO  TO  30 

CALCULATE  THE  GRCUF  AND  FC3ITI0N  INDEPENDENT  CONSTANT 


M 
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T4=T2/2. 

T5=T4*T  3 

T6=T2**? 

T7=2.*T6 

TB=2.*T7 

T9=T?**2 

T1Q=T9*2. 

Tll=T3/2. 

T16=H**2 


\ 


START  THE  CALCULATION  OF  THE  GROUP  FLUENCES. 


GROUPED 

4 J  ;  GR0UP=GR0UF+1 
PFWTMO  1-” 

IF (GROUP .GT .  NGRCUP) GO  TO  50 

CALCULATE  THE  GROUP  VIRGIN  FLUENCY  AND 
THE  GROUP  SCURCR  (S)  MATRIX. 


SCURCE=YIELO*SPECT(GROUD) 

^iLt,V3ERgJtiS5U2£E»PI*H»T16*X7NOB»X3^«^^IN.NUP,NHOR,l 

ADELX3,0,SIGS,S|jELR,X7SW,SIGTR,SIGR,SIGT,GR0UF,NGR0UP) 


SET  jo  T i J E  SUSMATRIC^S  DIA, UPPER,  ANO  LOWEP 
IN  PACKED  FORM. 


FIRST  CALCULATE  DCSITI0N  INDEPENDENT  CONSTANTS 

T1=T16*SIGR/D 
T 15=1 , /DEL0 
T17=Ti5*T15/2, 

T18=T15**2 
T19=T18*T16 
T'?G  =  ?.*T19 
T21  =  2.*T2‘' 

T22=Ti+T2. 

T12= 1 . / ( ?«  *0) 

T13=2.*Ti2*T3 
T14=T12*T2 
i  X?=X7MIN 

00  41  J=1,NUP 
C 

C  CALCULATE  THE  HEIGHT  CEPENCENT  CONSTANTS 

C 

71=EXP(XJ) 

Z2  =  EX°  T'l) 

Z7=Z1**2 

Z4=?2**2 

ZE=T5/Z1 

Z6=T7/Z4 


» OELX  J  ( 
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Z7=T it/Z3 
ZG  = T6/Z4 
Z9=T4/Z4 
ZlD=TS/?4 

IF <J.£Q.NUP.aND„L*EO. 1)3,2 

r-  i  .i  j  -  ^  o  /  1 

.  ■.-i  4  <7< 

Z13=Til/7o 
GO  TO  4 

3  Z14=T13/(Z3*Z2> 

Z15=T12/ (Z1+Z2) 

Z16=T12/<Z3*Z2> 

Z19=Z14+Z15-Z15 
Z17=T14/ <Zi*Z25 

4  IF(X3«GE«X3SW)G0  TO  1-0 

Z18  =  T 1+Z6+  Z7 
X1=X1MIN 

CALCULATE  01 A »  LOWER,  AND  UPPER 

00  14  I=1,NM0R 
X=Xl/75 
C  =  Z18 

IF(I. EG# 1) GO  TO  6 
Gl=-X 
G3=X 

C=C+T7* (Xl**2) 
9PAP=78+T6*(Xl**2)-79/Xl 
IF (I • cQ«  NHCR) GO  TC  5 
3?=X 
G4=-X 

3=Z8+T6* (Xl**2) +Z9/X1 
GO  TO  7 

5  8  =  0 
G2=C « 

G4=0. 

GO  TO  7 

6  Gi  =  3  • 

G2=°  • 

G3=0. 

G4=P . 

3PAR  =  i*  . 

3  =  710 
C=C+76 

7  IF<J.E0.NuP.AND.L.F0.1 >00  TO  11 

IFU.EG.DGO  7r  3 
APAR=Z11-Z12+713 
IF(J.EQ.NUF) 1C ,9 

8  AP  AR=  r; 

G1=0. 

G2=r'. 

9  A=Z11+Z12-Z13 
GO  TO  13 
A  =  % 

G3=C « 


10 
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G4=0. 

GO  TO  13  \ 

11  AsO. 

G1=G  « 

G2=0  • 

G3=u  • 

G4=0» 

APa:?  =  ?7 
C=C+Z13 

IFd.tO.DGC  TO  12 
9PAR-B9AR- (Xl/717) 

'f •  .’.fcfl«Nrfn«»)i;o  TO  IT 

12  d-bT  \  _ 

13  DIA(I,l>s90A» 

0IA(I,2)=-C 

0IA(I»3)=P' 

LOWER  (I»  1) =G1 
LOWER (1*2) =  ABAR 
LOWER  { I  *  3) =02 
UPPER (1,1) =G3 
UPPER (I»  2) *A 
UPPER (1,3) =  04 

14  Xi=Xl+DELXl 
0FLRL0W=QELX1*H*Z2 
GO  TO  120 

ICO  ■  IF<J.EQ.NUP,ANO.L.EO.i>GP  TO  113 
IF  ( J«  EG.NUF)  GO  TO  1U 
A=Zli*Z12-Z13 
IF(J.EQ.l) 11 J, 112 
11-  APAR=c. 

GO  TO  114 

111  A  =  Q. 

112  ABAo=Zll-Zl2+713 
GO  to  114 

113  ABA°=77 
A  =  C. 

114  Z18=Z7  +T  2? 

X1=0 

00  1„P  I=l,NMOR 
G  =  Z18 

IF(I.EO.l) GO  TO  132 
8BAR=T19- (T17/X1) 

IF(I.EO.NHCP>r-C  TC  131 
B=T19+(Tl7/Xl) 

GO  TO  103 
1 j  1  9  =  C  « 

GO  TO  1" 3 

132  C=C+T2j 
BRAR=-. 

B  =  T2l 

133  IF(J.EQ.NUP,ANC,L. EC. 1)1.4,105 

104  C=C+Z19 

105  OIA (I  *  1) =  9nAR 
0IA(I,2)s-C 
OIA(I-/3)=n 
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IF (X3. NE » X  3SW) GO  TO  1"6 
Kfl=i 

00  121  K=l,3 

121  L0W£R(I,K>=r. 

IF (I.Nt, 1) GO  TO  122 
KCOL ( I)  =  1 
LOWER  (T*i)  s/l^AR 
GO  TO  ir*  6 

122  IF(KA.£Q.?,0R.KA«£n,3> GO  TO  1D6 
K=Xi/DCLRLCW 

KK=K*i 

IF(KK.EQ.NHOR) KA  =  2 

IF(KK.GT.NH03)KA=3 

CHI=  (XI- ('<*C£LRLCW)  > /HARLOW 

LOWER ( 1,1) sA^AR-CHI^A^AR 

GO  TO  (123,124,124) KA 

123  LCWER(I,2)=CHI*A9AP 

124  KCOL < I ) =KK 
1 j6  X1=X1+0ELR 

c 

C  CALCULATE  THE  W  SUBHATPIX 

C 

120  00  15  I=1,NH09 

00  15  K=1,NH0R 

15  W  (K, I ) =u  .  0 

IF( J. Nc. 1) GO  TO  13 
W (1,1)  =  DIA  (1,2) 

W(1,2)=DIA  (1,3) 

W (NHOP , K1 ) =0 1  A ( NHOR , 1 ) 

W(NHOR,NHOS)=0IA(NHCR,2) 

K2=C 

00  17  1=2, K1 
00  16  K=1 , 3 

16  W(I,K+K2)=CIA  (I,K) 

17  K2=K2+ 1 
GO  TO  25 

18  IF(X3-X3SW)99,i3J  ,140 

99  00  19  1=1 , NHOR 

W(i,I)=»L0WER(l,2) ♦0(1,1) 

19  W (NHOD, I)=-LQW£R (NHCR,i)+G (Kl, I  ) -LOWER (NHO 9, 2) *0 (NHOR , I ) 
K2=0 

00  22  M-2, Ki 
00  21  1=1, NHOR 
SUM=n. 

•00  20  K=l,  3 

20  SUM  =  FUM  +  LOWlR(M,<)  *Q(’<  +  K2,I) 

21  W(M,I)=-5UM 

22  K2=K2+1 
GO  TO  131 

130  KA=i 

00  13^  H=1 ,NW0? 

IF(KA.E0.2.0?.KA.E0.3) GO  TO  16C 
KK=KCOL(M) 

IF (KK.t0#NH0o)KA=2 
IF(KK.GT.NHOR) KA=3 
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18..  00  175  I  =  1,NH0P 

suH=a.  ' 

GO  TO  (133,134,135)KA 

133  SUM=  LOWER(M»l)#Q(K<tr  ) +L0WER (M , 2) *Q (KK +1 »I) 

GO  TO  175 

174  SUM=  L0WES(M,1)*Q(KK,I) 

135  W(H,I)=-$UM 

GO  TO  131 

14)  00  141  M-l  ,NH03 

00  141  I  =  l,fW0R 
141  W(H»I)=-A3ARlfG(y»T) 

131  W (1,1)=0IA (1,2) *W  (1,1) 

W(i,2)=0IA (1,3) +W  (1 , 3) 

W (NHOn,Kl) =01 A (NFCP, 1 ) +W (NHOR, Kl) 

W(NH0P,NH0R)=CIA (NH0R,2)+W  (NHGF»NF0P) 

K2=C 

00  24  1=2, K1 
DO  2?  K=l, 3 

23  W(I,K+K2)=QIA(I,K)+W(I,K+K2) 

24  K2=K2+1 

CALCULATE  W  INVERSE 

25  CALL  MATRIX(10,NHCR,NHOR,C,W,6; ,Z7) 

CALCULATE  THE  Q  SU^MATRIX 

IF(J.EO.NUF) GO  TC  31 
IF(X7,GE,X3SW)GC  >'0  14? 

00  26  M=1,NH0° 

Q  (H ,  1 )  =  W(M,1)*UFPEP<1,2)+W(N,2)*L'FFER(2»1) 

26  0  ( M, NHOR)  =  W(M,K1)*UPFER(K1,3)+W()',NH0R)*UPFER(NHCR,2) 

PO  29  M=l,  NHO^* 

<2=0 

00  28  1=2, <1 
SUM=C. 

00  27  K=1 , 3 

27  SUM=SUM+W(M,K+K2) *UPPER(K+K2,4«K) 

Q  ( M , I )  =  SUM 

28  <2=<2+l 

29  CONTINUE 
GO  TO  144 

14?  00  14?  I=l,NHnR 

00  14?  <=iyNH0R 

143  Q  (I ,  K)  =?A*W(I,K) 

STORE  THE  0  SUPMATRIX  ON  CIS<. 

144  WPITE  (13)  ( (Q  ( I , K  )  ,I  =  l,NHCp) ,<  =  1,NH0R) 

CALCULATE  THE  G  SUB MATRIX 

31  IF(J.NE.i)00  TQ  34 
00  37  1= 1 , NHQR 
SL'M=0  * 


o  o  o 


GNE/PH/72- 8 


00  32  K=i,NHCP 

32  SUM=SUM+W(I,K)*S(K,J) 

33  G(I,J)=SU'< 

GO  TO  40 

34  JM=J-t 

IF(X3-X33W)i5u,151,i56 
15^  G<i,J)=LQWE9 <i,2)*G<l,JM) 

G (NHOR, J)= LOWER (NHOP,i)  *G(K1,JM>  ♦v.OWER(NHOR»2)*G(NHOR,J^) 
K2=C 

00  36  M  =  2  ,  K 1 
SUH=C'. 

DO  35  K=l,3 

35  SUM=SUM*L0WER<M,K)*Ge<+K2,JM) 

G (M, J)  =  SU‘* 

'36  K2=K2+1 

GO  TO  155 

151  KA  =  i 

DO  154  M=1 f NHOR 

IF(KA*EQ*2*0R*KA« £0*3)00  TO  131 
KK=KCOL ( M) 

IF(KK.EQ.NH0R)KA=2 
IF  <KK*GT.Nh0R)KA=3 
131  SUM=o. 

GO  TO  <152,153,154)KA 

152  SUM=LOWER < M, 1) *G ( KK, Ji) 4-LOWER (M, 2) *G<KK+1,JM) 

GO  TP  154 

153  3UM=L0WE9(M,1)*G(KK, JM> 

154  G(M,J)=SUM 
GO  TP  155 

156  00  157  M=1,NHPR 

157  G(M,J)=A9AR*G (M, JM) 

155  00  37  1=1, NHOR 

27  0IA(I,1)=S (I, J)-G (I,J) 

•00  30  1  =  1,  NHOR 
SUM=n, 

00  38  K=l, NHOR 

39  SUM  =  SUM+W(I,i<)*OIA  <K,1) 

39  G (I , J) =  SUM 

40  X3=Xo+0£LX3 

CALCULATE  FLUENCE  AT  EACH  POINT  IN  MESH 

JM=NUP+i 
00  48  J=i,NUF 
J J=JH- J 

IF( JJ.NE.NLP) GO  tq  42 
00  41  I=1,NH0P 

41  FLU ( I , J J)  =  G ( I , J J) 

GO  TO  46 

42  00  44  1=1, NHOR 
$UM=C. 

00  43  <=1,NHQP 

43  5L'M=SUM+  0  (  I ,  K 

44  FLU  ( I , J J ) =  SUM 
00  45  1=1, NHQP 


)  ’FLU  CK,  JJ4-1 ) 
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45  FLU(I,JJ>=G(I,JJ)-FLU<I, JJ) 

\ 

°EAD  A  G  SU9MATRIX  FRCM  DISK. 

46  IF( J.EG.l) GC  TO  48 
BACKSPACE  1C 
IF(J.E0.2)G0  TO  47 
9ACKSP ACE  1C 

47  READ  (1.)  < (C C I , K  ) , I=1,NH0R> ,K=l,NhOR> 

43  CONTINUE 

COMPUTE  THE  TCTAL  GROUP  FLUENCE. 

CALL  ACO(rLU,NUP,NHCRl 

PRINT  THE  DETAILED  GROUP  OUTPUT  IF  DE3IREC. 

IF (LOUT • EO • 0 ) GO  TO  49 

CALL  0UT<NUP,NH0R,X3MTNfXlMIN,CELX3,CELXifH,T|0ELRt X3SW,SIGT,GRCUF 
A) 

49  GO  TO  4 o'? 

CONVERT  THE  MULTIGROUF  FLUENCE  TC  A  TOTAL 
FLUENCE  OR  COSE  DEPENDING  UPCN  f'CCE  SPECIFIED. 

50  CALL  CONVFRT(NHCR,NUP,MCOE,NGRCUF,NNELT,NGAM) 

DETERMINE  IF  AIRCRAFT  SURVIVE. 

CALL  CHECK (NKCR, NUF, MODE, VUL, PCS, NUHPER,X3MIN,DELX3,X3SW,DELR , 
AX1MIN,9ELX1,H,ARRAY,L,RURST, ACFOS) 

DRAW  PLOTS  IF  DESIRED. 

IF (LMA ° a  EO , 0 ) GO  TO  51 

CALL  HaP(LFAPfNUF,rH0o,P0S,MCCF,VUL,HC5,DELR»>lMIN,DELXl,X3SW, 
AX?MIM,QELX3,L , ARRAY, NUMRER,H> 


51 

REWIND 

9 

°EWTNC 

1** 

RFWIMC 

11 

?oC 

RETURN 

END 
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Subroutine  SELECT 

\ 

This  subroutine  determines  what  weapon  output  spectrum 
will  be  used.  The  user  specifies  the  type  of  weapon  and 
the  source  of  the  spectrum.  He  can  load  his  own  spectrum 
or  use  one  of  the  unclassified  default  spectra  stored  in  the 
labeled  common  SPECTRA.  This  subroutine  is  called  by 
GAMNEUT 

Subroutine  SELECT  Glossary: 

GAMSO 

NEUTSO  A  number  that  determines  if  the  default 

spectra  or  a  user  supplied  spectra  is  to  be 
used 

NGROUP  The  number  of  neutron  plus  gamma  groups 

NNEUT  The  number  of  neutron  groups 

SPECFG  The  default  fission  gamma  spectrum 

SPECFN  The  default  fission  neutron  spectrum 

SPECTNG  The  default  thermonuclear  gamma  spectrum 

SPECTNN  The  default  thermonuclear  neutron  spectrum 

TYPE  The  weapon  type 
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SUB ROUTT  IE  SELECT (SF?CT, NNEUT,NGAM,NGRGUP) 

THIS  SUBROUTINE  Rr  ADS  IN  THE  WEAPON  TYPE, 

SOURCE  OF  WFAcon  SPECTRUM,  AND  SPECTRUM 
(GAMMA  AND  NEUTRON)  IF  SUPPLIED  BY  USER. 

WEAPON  TYPE  (TYCE> 

-  A  ONE  IS  A  MISSION  WEAPON 

-  A  TWO  IS  A  THERMONUCLEAR  WEAPON 

SOURCE  CF  WEAPON  ( GAMSO  FOR  GAMMAS, 

NEUTSO  FOR  NEUTRONS) 

-  A  BLANK  (ZERO)  INGICATES  TEE  CEFAULT 
SPFCTRUM  WILL  EE  UCED 

-  A  ONE  INDICATES  A  USER  SUPPLIED 
SPECTRUM  WILL  EE  USED 

.  MOTE..  IF  THE  DEFAULT  SPECTRUM  IS  NOT  USEC, 

THE  NUMBS3  OF  GROUPS  AND  ENFRGY  BANCS  MUST 
AGREE  WITH  THAT  S'JPPLIEO  AS  INPUT  FOR 
GROU°  CROSS  SECTIONS  AND  CCSE  CALCULATIONS. 

DIMENSION  5PECT<4:> 

INTEGER  GAMSC,  TYPE 

COMMON  GAMSO, TYPE, NEUTSO, I, M 

COMMON/ S RE CTRA/SP EC TNN  (22)  ,SF£CFN(22)  ,SPECFG(16)  ,SPECTNC-(13) 
REAO  1, TYPE, NEUTSC, GAMSO 

DETERMINE  SOURCE  OF  NEUTRON  SPECTRUM. 

IF(NEUTS0.E0.1)GC  TC  5 

DETERMINE  WEAPON  type, 

IF (TY°E. EO .2 ) GO  TO  4 

'  LOAD  DEFAULT  NEUTRON  FISSICN  SPECTRUM. 

00  3  1=1,2? 

SPECT(I)=SPECFN(I) 

GO  TO  7 

LOAD  OpF AULT  NEUT°CN  TN  SPSCTCUM. 

DO  5  1=1,22 
SPECT(I)=SPEOTNN(I) 

GO  TO  7 

LOAD  USER  SUPPLIED  NEUTRON  SCECTRUM 
READ  1j  ), (SPECT(I) ,I=1,NNEUT) 

DETERMINE  SOURCE  OF  GAMMA  SP^CTFUM. 


IF  (GAMSO. EQ.i) GC  TC  11 
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DETERMINE  WEAPON  TYPE. 
IF(TYF£.E0.2)GC  TC  2Q 

LOAD  DEFAULT  GAMMA  FISSION  SPECTRLM. 
00  21  1=1,18 

21  SP£CT(22+I)=S°ECFG(I> 

GO  TO  12 

LOAO  DEFAULT  GAMMA  TN  SPECTRUM. 

28  00  3  1=1*18 

8  SFECT<22*I)=SFECTNG(I> 

GO  TO  12 

.  LOAO  USER  SUPPLIED  GAMMA  SPECTRUM > 

11  M=NNEUT+1 

READ  lu:,(SPECT(I),I=M,NGROUP) 

1  FORMAT ( 1314) 

10"  F0RMATC1Y,1P7£H.u) 

12  RETURN 
ENO 
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Subroutine  MESH 

This  subroutine  calculates  the  number  of  vertical 
and  horizontal  mesh  points  and  the  spacing  between  mesh 
points.  The  maximum  number  of  horizontal  mesh  points  used 
is  40;  however,  the  number  can  be  changed.  The  card 

IF  (NIIOR .  GT.  40)NHOR=40 

limits  the  number  to  40.  If  any  other  number  of  mesh  points 
is  desired,  the  user  can  change  the  two  40*s  to  that  number. 
This  number  can  not  be  greater  than  60. 

The  maximum  number  of  vertical  mesh  points  is  70,  but 
this  number  can  also  be  changed.  The  card 

IF (M. GT . 70)M=70 

limits  the  number  to  70.  If  any  other  number  is  desired, 
the  user  can  change  the  two  70’ s  to  that  number.  This 
number  can  not  be  greater  than  120. 

If  the  number  of  mesh  points  are  increased,  the  running 
time  will  also  increase.  At  present,  this  option  has  not 
been  exercised,  so  no  estimate  of  the  increased  running  time 
can  be  made. 

Subroutine  MESH  Glossary: 

X1MAX  The  maximum  value  of  x* 

3 

X3MAX  The  maximum  value  of  x 

X3T0TAL  The  total  length  of  x^ 

See  the  GAMNEUT  glossary  for  the  remaining  terms. 
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SUBROUTINE  V-SH{H,STMTN,STMAX,H08,X1WIN,DELX1,NH0R,1,X3MIN» DELX2, 
A«JP,X3H0B> 

COMMON  Z,I,J,T1,K,m,B,N,X2TOTAL,X1MAX,X3MAX,TT,T 
C 

C  THIS  SUaRCUTINE  CALCULATES  the  number  of 

C  VERTICAL  AMO  HORIZONTAL  MESH  FOINTS  ANC  THE 

C  SPACING  BETWEEN  M-SH  POINTS.  IT  ALSO 

C  GFTERii'NES  THE  MINIMU''  AND  M  AXI MUM  FOR  CACH 

C  COORDINATE.  CHECKS  APE  HADE  TO  PREVENT  THE 

C  TOP  OF  THE  HE SH  BEING  GREATER  THAN  l-j“  KM 

C  AMO.  THE  BOTTOM  CF  THE  MESH  BEING  LESS 

C  T«EN  OP  EQUAL  TO  G  KH.  IN  ADDITION  A  m£SH 

C  POINT  IS  NOT  ALLOWED  TO  3E  ON  THR  HOB. 

C 

C  X1MIN  IS  THE  MINIMUM  XI  COORDINATE. 

C 

C  X1MAX  IS  THE  MAXIMUM  Xi  COORDINATE. 

C 

C  D^LXl  IS  THE  XI  INTERVAL. 

C 

C  NHQ3  IS  THE  NUMBER  OF  MESH  POTNTS  IN  THE 

C  XI  OIRECTTON. 

C 

C  L  IS  7ER0  UNLESS  THE  TOP  OF  THE  MESH  IS 

C  lu2  KM.  1^  THE  TOP  OF  THE  MESH  IS  ICO  KM 

C  l  IS  SET  TO  1. 

C 

C  X3MIN  IS  THE  M  INTHUM  X3  COORDINATE. 

C 

C  X3MAX  IS  THE  MAXIMUM  X3  COORDINATE. 

C 

C  DELX3  IS  THE  X?  INTERVAL. 

C 

C  NUP  IS  THE  NUMBER  OF  MESH  POINTS  IN  T»E 

C  X3  -DIRECTION, 

C 

c 

C  XlMIN,XlMAXfNHPC,  fi^D  DELX1  ARE  CALCULATED 

C  FIRST.  X1MAX  IS  EQUIVALENT  TO  1C  MEAN 

C  FREE  PATHS  CF  THE  GROUP  WITH  THE  LONGEST 

C  MEAN  FREE  PATH. 

C 

TrH'STMIN 

TT=H*STMAX 

X1MIN=1. 

X1MAX=1G ./T 
NHOP=o*TT»X1MAX 
IF(NH0k.Gt.‘*',)NHCF  =  4T 
DELX1=X1MAX/NHQR 
C 

C  Nr XT  XImxn  is  LCCATEO  ll  MEAN  FFEE  Paths 

C  DOWN  FROM  HCS.  WOH-VER,  7  IS  NCT  ALLOWED 

C  TO  BE  EOUAL  TC  Ca  LESS  THEN  0  KM. 

C 

Z  =HO° 
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fe 


i 


ij 


f 

f 

t!  j 


I 


II 


ft 


00  1  1=1, 10 
J=I 

T1=H*AL0G(1.h(£XB (Z/H)/T) ) 

Z=Z-T1 

IF (Z.GT.u) GO  TO  1 

Z=7+T1 

GO  TO  2 

1  CONTINUE 

2  X2MIN=AL0G<Z/U> 

C 

C  NeXT,  X3MAX  IS  LOCATEC  10  MEAN  FREE  PATHS 

C  UP  FROM  HpR.  HCWEVEP,  Z  IS  NCT  ALLQWEC 

C  TO  PE  GREATER  THEN  lw  *.  K-.  IN  AQCITIGN, 

C  A  CHECK  IS  MACE  F0=>  AN  INFINITE  MEAN  FREE 

C  PATH  IN  WHICH  CASr  Z  IS  SET  TC  1C'  KM. 

C  IF  7  IS  SET  TC  1j'  KM,  L  IS  SET  TC  1. 

C 

L  =  C 
Z=HOR 

HO  5  1=1,1 j 
K=I 

T1=FXP(Z/H)/T 

IF ( T1 • LT. 1 • ) GO  TO  4 

3  L  =  1 
Z=l.dE+G7 

’  GO  TO  6 

4  Z=7-H*AL0G(1.-T1) 

IFd.GE.i.  OE  +  17)  GO  TO  3 

5  CONTINUE 

6  X3MAX  =  AL0G  (7/H) 

C 

C  NEXT,  DELX3  ANO  NUd  ARE  CALCULATES.  A 

C  CHECK  IS  MAGE  TC  INSURE  THAT  THE  HC3  IS 

C  NOT  ON  A  MESH  POINT. 

C 

J=J  +  K 

M  =  3.0*J*ALCG(1.  +  (EXF(HOE/H)/T) ) /ALOG ( 1 .+  (EXP ( HC 5/H)  /TT)  ) 
IFtM.GT.  70)M=71 
X3T0TAL=X3^AX-X2MIN 
X3H0E=AL0G (HOC/H) 

T1=X3H03-X?MIN 

7  0ELX2=X3TCTAL/M 
Q=T l/nELX7 

N=B 

IF ( (8-N) .NE. 3 ) GO  TO  3 
M=M-1 
GO  TO  7 
3  NUP=M 

X3MIN=X3MIN+0ELX3 
IF  ( L  .  EO .  1)  GO  TO  <3 
NUP=NU?-1 
3  RETURN 

ENO 
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Subroutine  LOCATE 

This  subroutine  checks  the  aircraft  positions  against 
the  mesh  area  set  up  in  subroutine  MESH.  If  all  the  air¬ 
craft  are  outside  the  meshed  area,  NRET  will  remain  zero  and 
cause  subroutine  GAMNEUT  to  terminate  fluence  calculations. 
If  at  least  one  aircraft  is  inside  the  meshed  area,  NRET  is 
set  equal  to  one  and  fluence  calculations  will  continue. 

Subroutine  LOCATE  Glossary: 

A  The  scale  height  of  the  atmosphere  in  kilo¬ 

meters 

ACPOS  An  array  that  stores  aircraft  positions  in 

X,  Y,  Z  coordinates 

ARRAY  An  array  that  specifies  if  the  aircraft  is 

still  live 

BURST  The  coordinates  of  the  weapon  burst 

DELX1  The  x*  interval 

DELX3  The  x^  interval 

H  The  scale  height  of  the  atmosphere  in 

centimeters 

NHOR  The  number  of  horizontal  mesh  points 

NRET  A  number  that  indicates  if  any  aircraft  are 

in  the  meshed  area 

NUMBER  The  number  of  starting  aircraft 

NUP  The  number  of  vertical  mesh  points 

POS  An  array  that  stores  the  aircraft  positions 

in  r,z  coordinates 

A' 

R  The  radial  distance 

XI  The  value  of  x* 

X1MIN  The  minimum  value  of  x* 
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X3MIN 


ZMAX 


ZH  IN 


The  minimum  value  of 

The  altitude  of  the  upper  boundary  of  the 
meshed  area 

The  altitude  of  the  lower  boundary  of  the 
meshed  area 


105 


oooooooooooo 


— V  ■**» 


W'V- 


GNE/PH/72-8 


SUBROUTINE  LGCATE (NUMBER, ACPOS,8URS7,X1MIN,DELX1,NHOR,H,X?MIN,OELX 
AT, HUP, L,POS,NpET, ARRAY)  \ 

THIS  SUBROUTINE  LOCATES  THE  AIRCRAFT'  IN 
A  R,Z  GECNFTRY  NI7M  RESPECT  TO  T-»E  BURST. 

THE  AIRCRAFT  POSITIONS  ARE  CHECkEC  AGAINST  THE 
OUTER  90UN0ARIFS  OF  THE  mcsh  TO  SEE  IF  ANY 
AIRCRAFT  ARE  W I T H I f 4  THE  MESH.  CNEE  AN 
AIRCRAFT  IS  LCCAT £ J  INSIDE  THE  HESHEO  AREA 

this  subroutine  terminates  anc  control 

RETURNS  TO  GAMNE'JT.  IF  ALL  THE  AIRCRAFT 
ARE  OUTSIO-  THE  MESHEC  AREA,  KRET  REMAINS 
ZERO. 

DIMENSION  ACPOS (3,1.0) , BURST  CM  , PCS <2  ,100) 

INTEGER  ARRAYTIjO) 

COMMON  A,9,X1,ZHIN,ZMAX,R,T 
A=H*1..jE--»5 
NpET=f 

no  10  1=1 , NUM?E° 

IF<ARRAY(I).EQ.3)GO  T"  10 

°oS(2,I)=ACP0Sn,I) 

POS(i,I)=SCRT { ( (BURST (1) -ACPCS(1,I) )**2)  +( (BURST (2) -ACPCS (2, I ) ) **2 
A)) 

10  CONTINUE 

ZMIN=A*£XP (X3MIN-OELX?) 

9=-0ELX3 

IF(L.NE.1)S=0. 

ZMAX=A*EXP  (X3MIN+  (N'UP“OELX3)  *°) 

Xi=NHOR*QELXl*A 
DO  11  I  =  l,NiJMorp 
IFtAPRAY  (I)  .EO.",)GO  Tn  11 

IF (PCS (2, I) .GF.ZHAX.OR.POS  (2, I) . LE.ZNIN)  12,13 
13  NRET=1 

GO  TO  2n 

12  R=X1*EXP<  PCS (2 , I ) /A ) 

IF(POS(l,I)  «GE.R.CR.Pr'S(l»I)  .GE .  4C  0  .  )  11, 13 

11  CONTINUE 
2!  RCTUPN 

END 
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Subroutine  VIRGIN 

\ 

This  subroutine  has  three  functions.  First,  the  group 

virgin  fluence  is  calculated  at  each  point  in  the  mesh. 

Second,  a  check  is  made  to  determine  if  the  horizontal  radius 

exceeds  200  kilometers.  If  the  radius  is  greater  than  200 

3  3 

km  at  altitude  x  ,  the  value  of  x  is  stored  in  X3SW,  and 
the  coordinate  system  is  switched  from  the  nonorthogonal 
system  to  the  orthogonal  system.  The  third  function  is  the 
calculation  of  the  S  matrix  for  the  matrix  equation  AF  =  S. 
The  virgin  fluence  is  stored  on  a  disk  file  for  later  use. 

Subroutine  VIRGIN  Glossary: 

D  The  diffusion  coefficient 


DELR 

The  actual  radial  spacing 

DELX1 

The  x1  interval 

DELX3 

3 

The  x  interval 

DELZ 

The  distance  between  two  vertical 

mesh  points 

DELZSQ 

DELZ  squared 

F 

Group  total  fluence 

GROUP 

A  counter  relating  for  what  group 
fluence  is  being  calculated 

the  virgin 

RH0SQ 

The  distance  between  the  burst  and 
squared 

a  point. 

VIRG 

The  group  virgin  fluence 

XSECT 

The  group  cross  sections 
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SUPRCUTINF  V IRGIN (SOURCE, FI, H,TC,X3HCE,X3MIN, XI  MIN, NUP,NHOR,T,OSLX| 
Al.,0FLX3,T,SiGS,S,CELR,X3SW,Sir-TR,SIGR,SIGT,GRCUF,NGR0l'3)  | 

DIMENSION  S(B',123),F(6j,120)  j 

COMMON  Ti,T2,T3,X3,J,T4,T5,TS,OELZ,CELZSQ,T9,  Xi , I, RHOSC , VIRG ( 6 u ,12 j 
AO)  ,RA,K,T3,XSECT (43),L,LL,M  3 

INTEGER  6°  CUP  3 

EQUIVALENCE  <VI=G,F)  I 

THIS  SUBROUTINE  CALCULATES  THE  VIRGIN  1 

FLUENC6  AT  EACH  MFSH  POINT  AND  STORES  THE  1 

RESULTING  MATRIX  CN  TAPE.  THE  S  MATRIX  | 

FOR  THE  MATRIX  EQUATION  AF=S  IS  CALCULATED 

AND  RETURNED  TO  T*E  MAIN  PROGRAM.  J 

M=NGR0UP+3 

RtAO  (9)  (XSECT(I) ,I=i,M) 

SIGTR=XSECT( t) 

SIGR=XS£CT<2) 

SIGT=XS£CT (3) 

STGS=XSECT (4) 

D=i • / ( 3 • *S IGTR) 

T=H*STGT 

IF (SPU°CE. EC. C«  J) 16,17 

16  00  18  J= 1 , NUP 
DO  18  1=1 , NHCP 

18  V-IRG  ( I ,  J)  =  T  • 

GO  TO  19 

17  K  =  1 

T8  =  ^0URCE/  (4,*PI) 

T1=T8/TQ 
T2=EXP ( X3HCP) 

T3=l . /EXD ( T2) 

X3=X3MIN 

00  4  J=1,MUP  1 

T4=EXP(X3> 

TE=EXD (T4)  j 

T6=l • /T5  I 

T9=T6-T3 
0ELZ=T4-T? 

IF(G°CUP.ED.l) 7, 3 
«  IF(X3.GE.X2SW)2,R 

7  IF(K»  EQ. 2) GO  TO  2 

9  0ELZSQ  =  DEL  Zt*2 

X1=X1MIN 
DO  1  1=1, N^O3 
RH0S0  =  DEIZSCM(X1*T5) **2) 

VIRG(I,J)=(T1/RHCSC)*EXP( (T»SQRT(RHCSC)/DELZ) *T9) 

1  X1  =  X1  +  DELX  1 

IF(GP0UP.NE.1)G0  TO  4 
X3SW=1.  .Etc 
RA=H*X1*T5 

IF(RA.LT.2. JE+07) GO  TO  4 
DELR=RA/NHCR 

K-Z 

X3SW =X3 


\ 

i 
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GO  TO  4 

2  OELZ=H*DElZ 
OELZSQ=DSLZ**2 
RA=fl- 

00  3  I=i,NHQR 
RH0S0=0ELTSQ  +  (C?A*¥2) 

VIRG  (I ,  J  )■-  (T3/RFCSQ)*  EXP(  (T*SGRT (RHCSC) /OELZ)  *T9) 

3  RA=RA  +  0EI_R 

4  X3=X7+0£LX  3 

19  T1=(-TQ*3TGS) /O 

WRITE-  (il)  ( (VIRG  (I,J) ,I  =  i,NH0R) ,J  =  1,NUP) 

REWIND  11 
00  5  J=i ,NUP 
GO  5  1=1 *NHOR 

5  S(I,J)=T1*VIRG<I,J> 

I r(GP0UPtrQ»l)G9  TO  Z. 

K=GR0UP+4 
J=GR0UP-1 
00  ICi  1  =  1, J 
K=K~1 

IF(XSFCTC<)  .EO.  3)  11,12 

11  PEAD  (11) 

go  to  n, 

12  Tl=(-Td»XSECT(K) )/0 

READ  (11)  ( (F(L,LL) ,1=1, NHOR) » LL=1 » NUP) 

■  00  13  1.1  =  1, NUP 
00  13  L=l, NHOR 

13  S(l,LL)=S(l,ll) +Tl*F(l,ll) 

10  CONTINUE 

2u  RETURN 

END 
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Subroutine  ADD 

The  function  of  this  subroutine  is  to  add  the  group 
virgin  fluence  to  the  group  scattered  fluence  to  produce 
the  group  total  fluence.  The  group  virgin  fluence  is  read 
from  a  disk  file  and  the  group  total  fluence  is  written  on 
the  same  record  of  this  disk  file. 

Subroutine  ADD  Glossary: 

FLU  At  first,  the  group  scattered  fluence  and 

later,  the  group  total  fluence 

NHOR  The  number  of  horizontal  mesh  points 

NUP  The  number  of  vertical  mesh  points 

VIRG  The  group  virgin  fluence 
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SUBROUTINE  AQD( FLU,NUP,NHCR) 

\ 

THIS  SUP°OUTTNE  ADOS  THE  VIRGIN  ANC 
SCATTERED  FLUENCE  FOR  EACH  GROUP  TC  GIVE 
THE  TOTAL  GRCUP  FLUENCE. 

DIMENSION  FLU(6 3,12c) 

COMMON  I,J,VIRG(6-,12-) 

READ  (11)  ((VIRG(IyJ) ,I=1,NHCR) ,J=i,NLF) 

DO  1  J=1,NUP 
DO  i-  I  =  l,NHOR 
1  FLU (I, J) =FLU (I, J) ♦VIRG(IfJ) 

BACKSPACE  11 

WRITE  (11)  ( (FLU(I,J) ,I=1,NHCR) ,J-1,NLP) 

END 
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Subroutine  OUT 

This  subroutine  prints  a  detailed  output.  A  fluence 
for  every  group  is  printed  for  every  mesh  point.  This 
subroutine  is  called  from  GAMNEUT  only  iv  LOUT  is  equal  to 
one.  Since  the  output  will  be  several  hundred  pages,  the 
use  of  this,  subroutine  is  not  recommended  unless  the  user 
wants  to  examine  the  group  fluences. 

Subroutine  OUT  Glossary: 

A  The  scale  height  of  the  atmosphere  in  km. 

DELR  The  actual  radial  mesh  interval 

DELX1  The  x*  mesh  interval 

3 

DELX3  The  x  mesh  interval 

FLU  ’  The  group  total  fluence 

GROUP  The  group  number 

R  The  actual  radial  distance 

X3SW  The  value  of  x3  when  the  change  in  coordinate 

systems  occ»  r 

Z  The  altitude 
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SUBROUTINE  OUT ( NUF, NHOR, X3MIN,X1MIN , CELX3, OEL XI ,H, OELR,X3SW,GROI 
COMMON  FLU  (6 12?) ,  A,  X3,  I,  J,  9,  Z,X1,R  (60 ,C,K,M, JJ,G,RA 
INTEGER  SpOUP 
C 

C  THIS  SUTRPUT INE  PRINTS  THE  OETAILEC  OUTPUT, 

C  A  FLUENCE  IS  PRINTED  FOR  EVERY  MESH  POINT 

.  C  IN  EVERY  GpOUF. 

C 

A=H*l.nE-"5 
X3=X3MIN 
BACKSPACE  11 

RpAO  (11)  ((FLU  (I, J) ,I=1,NHCR) ,J=i,NUF> 

PRINT  5j?, GROUP 
00  S  J=i,NUP 
B=EXP(X3> 

G=FXP(B> 

Z=A*R 

IF(X3.GE«X3SW)G0  TO  2 

X1=X1MIN 

C=A»G 

00  1  1=1, NHOR 
R(I)=C*X1 

1  X1=X1+DELX1 
GO  TO  4 

2  RA=C « 

00  3  1=1, NHOR 
R(I)=RA*l.TE-?5 

3  RA=RA+OELR 

4  K=1 
M=8 

00  5  1=1,30 
DPINT  5?i,Z 

PRINT  5C  3 » (F ( J J) , JJ=K,M) 

PRINT  5?5,  (FLU( JJ,J) , JJ=K,M) 

K=K+8 

IF(K.GT.NHCR)GC  TC  6 
M=M+8 

IF(M,GT.NHOR)M=NHCR 
CONTTNU*7 
Xc=X3+0cLX3 

)3  FORMAT (///,3X»*THE  GROUP  IS  *,T4,//) 

01  FORMAT (//,3X, 4 HEIGHT  IS  *,pi:.5,4  KM • 4 ) 
n  3  FORMAT (3X,4RA0IUS  (KM) 4,4X,8P1H«3) 

,5  FCpMAT (3X, ’GROUP  RLUENCE*,lp8E14»4) 

END 
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Subroutine  CONVERT 

This  subroutine  converts  the  multigroup  fluences  to 

the  units  specified  by  the  MODE  parameter.  At  the  present 

time  MODE  can  have  a  value  from  one  to  eight.  The  meanings 

of  these  values  are  described  in  the  subroutine  listing. 

Subroutine  CONVERT  Glossary: 

F  The  group  total  fluence 

FLU  The  result  of  the  conversion,  either  total 

fluence  or  dose 

MODE  The  units  of  FLU 

SILGAM  The  multigroup  gamma  response  functions  for 

silicon  dose  in  rads  . 

SILNEUT  The  multigroup  neutron  response  functions 

for  silicon  dose  in  rads 

TISGAM  The  multigroup  gamma  response  functions  for 

tissue  dose  in  rads 

TISNEUT  The  multigroup  neutron  response  functions 

for  tissue  dose  in  rads 
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SUBROUTINE  CONVERT (NrlQR,NUF,MOGE,NGRCLF,NNEUT  ,NGAM) 

THIS  SUBROUTINE  CONVERTS  THE  NULTIGROUF 
FLUENCE  TO  THE  MQCE  SFECIFIEO  IK  THE 
VULNERABILITY. 

MOD£=l  TOTAL  NEUTRON  FLUENCE  (N/CM2) 

M0DE=2  TOTAL  0AmMA  FLUENCE  (G/CN2) 

M032=3  NEUTRON  TISSUE  DOSE  (RACS) 

MOOE=<<  GAMMA  TISSUE  DOSE  (RADS) 

H00E=5  NEUTRON  +  GAMMA  TISSUE  COSE  (RAOS) 

MGQE=6  NEUTRON  SILICON  CCSE  (FACS) 

'•*OOE=7  GAMMA  SILICON  DOSE  (RACS) 

MOOE=!  NEUTRON  ♦  GAMMA  TISSUE  CCSE  (RAOS) 

DIMENSION  MODE ( 2) 

COMMON  SILK£UT(22) ,SILGAM(18)  ,TISNELT(22)  ,TISGAM(18)  ,  K, F  (6.>,  120 
AGROUP.I, J,FLU(6J, 120) .LIMIT 
INTEGER  GROUP 
REMIND  1J 
REWIND  11 

READ  (S)  (SILNEUT  (I) , I=1,NNEUT) 

READ  (9)  (SILGAM(I) ,I=1,NCAM) 

READ  (9)  ( TISNEUT (I)»I=1»NNEUT) 

REAO  (9)  (TISGAM(I)  ,I  =  1,NC-AM> 

'  K=MOOE(l) 

3  DO  4  J=1  * NUP 
00  4  1  =  1  * NHOR 

4  FLU(I,J)='.r- 

10  GO  TO  (1,2, 8, 12, 8, 19,19, 15)K 

1  GR0UP=1 
LIMIT=NNEUT 
GO  TO  5 

2  0R0UP=1 
LIMIT=NGAM 

5  IF(GPOU°.GT. LIMIT) Gn  T0  7 

READ  (11)  ((F(I,J),I=1,NHCR) .J-l.NUF) 

00  6  J=1,NUP 
00  6  1=1, NHOR 

6  FLU(I,J)=FLU(I,J)+F(I, J) 

GR0UP=GR0UF+1 

GO  TO  5 

7  WRITE  (IS)  <(FLU(I,J), 1=1, NHOR) ,J=1,NUF) 

GO  TP  51 

9  GR0UP=1 

L IMI T=NMEUT 

9  IF(GROUP.GT. LIMIT) GO  TO  11 

RFAO  (11)  ((F  (I,J),I=1,NHC°) ,J  =  1,NUF) 

00  3G  J= 1 , NUP 
DO  3C  1  =  1*  NWGD 

30  FLU ( I , J) =TISNCUT ( GPGUP ) *  F  (I*J) +  FLU  ( I »J) 

GR0UR=G?0UP+1 
GO  TO  9 

11  IF(K.t0.7) GG  TO  7 

12  GP0UP=1 
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LTMIT=MGAM 

14  IF(GR0UP.GT.LIMIT)C-0  TO  7 

PEAO  Cll)  ((F (I,J),I=i,NHCP) ,J=i,NUF) 
00  13  J=l» NIIP 
00  13  1=1, MHO*5 

13  FLU(I,J)=TISGAM(GROUP)*F(I,J)+FLU(T,J) 
GROUP  =GR0UP*1 
GO  TO  14 

15  GP0UP=1 
LIMIT=NNEUT 

16  IF(GR0UP.GT.LIMIT)C-0  TO  18 

READ  (11)  ((F(I,J),I=1,NHCP),J=1,NUP) 

00  17  J=1,NUP 
00  17  1=1, NHOR 

17  FLU(I,  J)  =SILNE'JT  ( GROUP )*F  (I,  J)  +FLU(I, J) 
GR0UP=GR0UP+1 

GO  TO  16 

18  IF (K.EQ*  6) GO  TC  7 

19  GR0UP=1 
LIMTT=NGAM 

20  IF(GROUP.GT.LIMIT)GO  TO  7 

READ  (11)  ((F(T,J), 1=1, NHOR) , J=1,MUP) 

00  21  J=l, NUP 
00  21  1=1, NHCR 

21  FLU(I,J)=SILGAM (GROUP) *F (I ,J) +FLU (I , J) 

•  •  GR0UP=GR0UP+1 

GO  TO  20 

5C  GO  T0( 51 *60, 51, 60, 60, 51, GO, 60 >K 
51  K=HOOF (2) 

GO  TO  3 

60  RFWINO  1U 
REWIHO  ii 

RETURN 

ENO 
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Subroutine  CHECK 

\ 

This  subroutine  places  the  aircraft  in  the  meshed  area 
and  calculates  the  neutron  and  gamma  fluence  of  dose  the 
aircraft  experienced.  This  calculation  is  compared  to  the 
aircraft  vulnerability  to  determine  if  the  aircraft  survived. 
Subroutine  NOTICE  is  then  called  to  give  the  printed  output 
concerning  the  aircraft. 

Subroutine  CHECK  Glossary: 

The  aircraft  positions  in  X,  Y,  Z  coordinates 
The  burst  coordinates 

The  gamma  level  (fluence  or  dose)  experienced 
by  the  aircraft 

The  neutron  level  (fluence  or  dose) 
experienced  by  the  aircraft 

The  aircraft  gamma  and  neutron  vulnerability 


ACPOS 

BURST 

GLEVEL 

NLEVEL 

VUL 
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SU3R0UTINE  (  HECK  CNHCR,NUP,MOCE,VUL*  PCS, NUM9ER,X3MIN,0ELX3,X3SV«, 
40ELR,XlMIN,0ELXl,M,/looaY,Lt3UPST,dCFCS) 


THIS  SUBROUTINE  LOCATES  THE  AIRCRAFT  AND 
DETERMINES  IF  GAMMAS  CR  NEUTRONS  HAVE 
KILLED  ">»E  aircraft 


REAL  MLEVEL 
INTEGER  ARPAfdjO) 

DIMENSION  MOOE(R)  ,VUL(2)  ,P0S(2,10o)  ,Bl!RS”(3),ACPOSC3,li?0) 

COMMON  ILCE)  ,FG  (6C  ,  12?)  ,  K,ZNIN,ZHAX,  N,ZA  ,S ,NL  , NH ,3, DEL ,  Zl, Z2, 

AXlMftX,RMAX,31KAX,R2MAX,NP,D,01,NF.l,FNl,FN2,FGl,FG2,NLEVEL,GLEVEL, 
RW,H1,W3,«4,FM(6j,120) 

OEL=  DELPHI  • CE-35 
S=H*i  ,  3E-’-  5 
XlMflX=NHOR*OELXl 
RNAX=NHOR*CEu 
K=2 

IF(MCP£(1> .EQ.5.CR.M00t(l) .EQ.g)K=l 
READ  Cl.)  ((FN(I,J) ,I=1?NH0°) ,J=1,NUP> 

IF(K.EQ.l>GO  TO  1 

READ  (1?)  <(FG(t,J) ,I=t,NHCR> ,J=1,NLP) 

1  ZMIH=X3MIN-0ELX3 
7MAX=(NUP*0ELX3)+7MIN 
IF(L.dQ.l) ZMAX=2MAX+0ELX3 
Nil 

2  IF<N.GT.NUMPED)GC  TC  10 
IF(APRAY<N).FC.T)GO  TO  9 
A=P0S(2,N) 

ZA=ALCG(A/S) 

IF(ZA.GE.ZMAX.OR.ZA.Lr.ZMIN)GC  TC  8 
NL=(ZA-7MIN) /OELX3 
NH=NL+1 
9=P0S ( 1, N) 

Zi=7MIN+(ML*CELX3> 

Z2=Zl+OELX3 
IF(Z1.GE.X3SW)GC  TO  3 
R1MAX=S»X1MAX*EXC (EXP (71) ) 

IF(Z2.Gc.X3SW)GO  TO  4 
R2MAX  =  S*Xlf*AX*EXF  (EXP  (Z2) ) 

GO  TO  5 

3  RlMAX=pMAX 

4  R2MAX=RMAX 

5  IF(plMAX*LE»P)FNi=0. 

IF(p2MAX,LE»T)GC  TO  8 
IF  (FNi.EO.  J.  )  r-Q  TO  6 
0=P1MAX/N«CR 

NP= (o/O) +1 
IF(NL,E0.C)FN1=3. 

IF(FNl.EO.T) GO  TO  6 
IF  {NR.EG*N^Or,)Vll=  •)* 

IF(K1.E0,.,)G0  TC  2C 
Hl=FN(NR+l*NL) 

21  WsFN(NR«NL) 

FNi  =  W*((9-  (NP*D> )*(W1-W)/D) 
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6  Q1=R2MAX/NH09 
NRl=(P/ni) *1 
IF(NRi.EC.NH0R)W4=C. 

IF(M4.cQ.’.)GG  TC  21 
W4=FN(MR<1,NH) 

21  W3=FN(tRftll>) 

FM2=Ko+C  (3-(?JRl*01)  )*  (W4-H3)/ni) 

NLEVFL=FNi+(ZA-21) *<FN2-FNi> /0ELX3 
IFCNLEVEL-VULCi) ) 23,22,22 

22  ARRAY  C M)  -?. 

GO  to  24 

23  ARRAY C  N) =4 

24  IFCARoay (N).cO.2.QR.<.E0.1)G0  TO  3 
IFCFHi.EQ. i.)FGl=i. 

IFCFG1.E0. J.)GO  TO  7 
IFCW1.E0.- .)G0  TC  25 
W1=FG (NR+1 jNL) 

25  W=FG<nR,NU 

FGi=W+< (9-(MR*0) )*(W1-W)/C> 

7  IF  CW4. £9 • ' • ) GO  TC  25 
W4=FG(NRl+i,NH> 

26  W3=FG (NPl  * NH) 

FG2=W3f( (B-(NP1»01) )* (W4-N3) /Dl) 

GLEVFL=FGl+ ( (ZA-Zl) *(FG2-FG1> /CELX3) 

IF  (GLEVEL-\/UL(2)  )2P,27,?7 

27  ‘  A PRAY (N) =3 

GO  TO  8 

28  ARRAY (N) =4 

8  CALL  NOTIC£(ARRAY,GLSVEL,  NLEVrL,  BURST  ,VUL,N,MCOF.,AC<=OS,K) 

9  N=N+1 

Wl=l. 

GC  TO  2 

10  RETURM 
END 
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Subroutine  NOTICE 

This  subroutine  is  called  by  subroutine  CHECK  and  its 
only  function  is  to  print  output  detailing  the  neutron  and 
gamma  levels  experienced  by  the  aircraft.  The  printout  also 
states  if  the  aircraft  survived  the  neutron  and  gamma  levels. 

Subroutine  NOTICE  Glossary: 

The  aircraft  positions 
The  burst  coordinates 

The  gamma  level  experienced  by  the  aircraft 
The  units  of  GLEVEL  and  NLEVEL 
The  neutron  level  experienced  by  the  aircraft 
The  aircraft  neutron  and  gamma  vulnerability 


ACPOS 

BURST 

GLEVEL 

MODE 

NLEVEL 

VUL 
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SUBROUTINE  NOTICE  (ARRAY, GLEVEL,NV-EVFL, BURST, VIL,N, MODE,  ACPOS, 

THIS  SUBROUTINE  HILL  PRINT  CUT  THE  RESULTS 
-  OP  THE  GAMMA  AND  NEUTRON  FNVT.ROMENT  ON  THE 
AIRCRAFT.  IF  THE  AIRCRAFT  HAS  BEEN  KILLED 
IT  WILL  STATE  AIRCRAFT  KILLED  ANC  GIVE 
BURST  AND  AIRCPAFT  LOCATION,  VULNERABILITY 
LEVEL,  ANC  ACTUAL  P-AMMA  OR  N.-U'TRCN  LEVEL. 

IF  THS  AIRCRAFT  SURVIVES,  PUT  DIC  FXPERIENCE 
SCHE  GAMMA  ANC  NEUTRON  LEVEL,  THE  PRINTOUT 
•HILL  STATE  THE  AIRCRAFT  SURVIVEC,  GIVE 
AIRCRAFT  AND  CURST  LOCATION ,  ANC  GIVE  ACTUAL 
GAMMA  OR  NEUTRON  LEVEL.  IF  THE  AIRCRAFT 
SURVIVED  AND  WAS  OUTSIGE  THE  MESH  AREA 
(THAT  IS,  THE  AIRCRAFT  DIC  NOT  EXPERIENCE 
ANY  GAMMA  OR  NEUTRON  LEVEL),  THE  FRINTGUT 
WILL  STATE  THE  AIRCRAFT  SUPVIVEC,  GIVE  THE 
.  AIRCRAFT  AND  BUPST  LOCATION,  ANC  STATE  THAT 
THE  EFFECTIVE  GAMMA  AND  NEUTRON  LEVEL  WAS 
ZERO. 

REAL  MLFVFL 
INTEGER  ARRAY (1j3) 

DIMENSION  BURST  (3)  ,V'JL  (2)  ,MO0E  (2)  ,ACPCS  (3,1,0  ) 

COMMON  I,M,L,J,KK 
M=MODF (1) 

L=MODE (2) 

J=ARRAY(N) 

KK=i 

PRINT  1,N, (ACPOS (I,N) ,1-1,3) 

GO  TO  (1., 11,12,13) J 

10  PPINT  2, (GUFST(I) ,I“1,3) 

PRINT  3 

GO  TO  53 

11  GO  TO  ( 14, 15) K 

14  PRINT  5, (BLPST(I) ,1=1,3) 

32  IF (M.EQ. 5) GO  TO  16 

17  PRINT  6 , VUL < 1 ) , NLEVFL 

GO  TO  53 

16  PRINT  7,VUL(l),Nl£VEL 

GO  TO  50  - 

15  PRINT  4, (BURST(I) ,1=1,3) 

3 j  IF(M.£0.3)GC  TO  17 

IF(M.FQ.6)G0  TO  15 
23  POINT  '3  ,  VUL  ( 1 )  *NL  EVEL 
GO  TO  53 

12  PRINT  51, (BURST (I) ,1=1,3) 

31  IF (L . EQ*  2 ) GO  TO  19 

IF(L.F0.4)G0  TO  22 
PRINT  6,  V'Jl  (2) ,  0LEVrL 
GO  TO  5j 

19  P°INT  9,VUL(2) ,GLEVEL 

GO  TO  5' 

22  PRINT  7 , VUL ( ? ) , GL EVEL 

GO  TO  5T 
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1?  GO  <23,2 1)K 

2d  PPTVT  5c,  *-URSTCI>tI=l,3> 

GO  to  32 

21  potnt  ?t (?LF?T(Ii f-  -.3> 

PRINT  53 
GO  TO  3“ 

53  IF(J.NE.4>G0  TO  1-0 
IFtKK.EO.P'GC  TO  lL" 

KK=2 

•VvlNT  54 
GO  TO  31 

1  FORMAT (////, 3X, ‘AIRCRAFT  NUMBER  *,I4,*  LOCATEC  AT*,F3.3,*  KM.,*, 
AF8.3,*  KM.,*,F3.?,*  <M.*) 

2  FORMAT (4X,*AIRCRAFT  SURVIVED  GAMMAS  ANC  NEUTRONS  FROM  BLRST  AT*, 

are. 3,*  km.,*,f8.?,*  km.,*, Fa.?,*  km.*) 

3  FORMAT (4X, ’EFFECTIVE  GAMMA  AND  NEUTRON  LEVrL , WAS  ZERO*) 

4  FORMAT (4X»*KILlEC  PY  NEUTRONS  FROM  8LRST  AT*,F8.3,*  KN.,*,F8.3,* 
AM.,*,F5.3f*  KM.*) 

5  FORMAT (4X, ’KILLED  EY  r'0M3INEC  GAMMA  +  NEUTRON  DOSE  FROM  BURST  AT* 
AF8.3,*  KM.,*,F8.3,*  KM.,*,F8S3,*  KM.*) 

5  F0RMAT(4X, ’VULNERABILITY  LEVEL  WAS* , 1FE 11. 4, *  RADS  (SILICON  DOSE) 
A  ACTUAL  LEVEL  WAS*,lPEii. 4,*  RADS  (SILICON  DOSE)*) 

7  F0RMAT(4X,*VULNERAFILITY  LEVEL  WAS*, 1cE11.4,*  PADS  (TISSUE  OCSE), 
AACTUAL  LEVEL  WAS* , 1PE1 1. 4, *  RADS  (TISSUE  COSE)*) 

8  FORMAT (4X, *VULNERAFILITY  LEVEL  WAS*, 1FE11.4,*  NEUTR0NS/CM2,  ACTLA 
A  LEVFL  WAS*,1CE11.4,*  NEUTR0NS/CN2* ) 

9  FCRMAT(4X,*VULNERARILITY  LEVEL  WAS*,1FE11.4,*  GAMMAS/CM2,  ACTUAL 
AEVEL  W AS* , 1PE11 , 4  ,*  GAMMAS/C*?*) 

51  F0RMAT(4X,*KILLEC  EY  GAMMAS  FROM  °URST  AT*,F8.3,*  <M .  , * , F 3 . 3 , *  KM 

A,*,Fg.3,*  KM.*) 

52  FORMAT (ux, ’AIRCRAFT  SURVIVED  COMBINED  GAM^A  +  NEUTRON  DCSE  FRCM  ? 
ARST  AT*,FS.3,*  Kw. CS. 3,*  KM.,*,F8.3,*  KM.*) 

53  FORMAT (4X, ’NEUTRONS*) 

54  FORMAT (4X, ’GAMMAS*) 

IB  return 

END 
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Subroutine  MAP 

This  subroutine  controls  the  plot  output.  It  calls 
subroutine  CONTOUR  to  calculate  the  isofluence  or  isodose 
lines  and  subroutine  SETUP  to  set  up  the  plots.  The  actual 
plotting  of  the  lines  is  done  in  this  subroutine. 

This  subroutine,  as  written,  is  designed  for  use  on  an 
on-line  plotter.  With  the  addition  of  two  cards,  however, 
this  subroutine,  and  therefore  the  entire  code,  can  be  used 
at  the  computer  center  where  only  an  off-line  plotter  is 
available.  The  two  cards  are: 

CALL  PLOTS (WORKA, 1024,7) 

CALL  PLOTE 

The  first  card  should  be  the  first  instruction  in  this  sub¬ 
routine  and  the  second  card  should  be  the  last  instruction. 

Using  this  option  will  require  different  control  cards 
since  a  magnetic  tape  is  necessary.  Therefore,  the  user 
should  check  the  local  instructions  before  this  option  is 
exercised.  Also,  the  first  card  of  the  MAIN  program  must 
be  changed.  The  file  called  PLOT  should  be  replaced  with 
a  file  called  TAPE7.  Tape  7  is  then  the  magnetic  tape 
required  for  the  off-line  plotter. 

Subroutine  MAP  Glossary: 

ALT  The  height  of  the  burst  in  km 

DATA  The  isofluence  or  isodose  data  calculated  by 

CONTOUR.  This  array  stores  the  r,z  coordi¬ 
nates  for  the  isofluence  or  isodose  points 
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SCALER 
SCALE2 
WORK  A 


The  x  axis  scale  factor  for  the  plot  routines 

The  y  axis  scale  factor  for  the  plot  routines 

An  array  that  is  unused  in  this  program.  If 
the  program  is  converted  for  use  on  an  off 
line  plotter,  this  array  will  be  the  work 
area  needed  for  the  plot  calculations 
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SU9R0UTIN-  -aP(LMAP,NUP,NHOR,PeS,MOCE,VUL»HOe,OELR,XiMIN,DELXl, 
4X3SW,X3MI.N,PELX3,L,A.RRAY,NAC  ,H) 

THIS  SUnRCUTINE  PLOTS  THE  GAMMA  AND  NEUTRON 
EMVIPOrENT.  THE  T YPE  OF  PLOT  IS  DETERMINER 
°Y  MODE  AND  LMAP. 

DIMENSION  M0DE(2) , VUL ( 2) , FOS (2 , 10 " ) 

INTEGER  APRAY(ljO) 

COMMON  aLClE)  ,0ATAC»1 '  ,16)  ,M(B)  ,p(6l  ,120)  ,RAD  , SCALER,  SC ALEZ,  K  ,M , 
ANPAS?,Kk,LIwIT,  JjLCC'JNT,  I ,NPD, ALT,X ,Y,6,N9,8A  ,NSTOP,NCCUNT,HI  ,C  ,D 
AJP,WORKA(102A) 

DETERMINE  HHAT  PLOTS  HILL  RE  CONE. 

REWIND  n 
ALT=HG=>*1.  IE-05 
MFASS=1 
N?TOP=2 
K=MOO£ (NPASS) 

IF(K.£0.5.CR.K.EQ.6>NST0P=i 

3-H*l < ji-. 5 

IF(X3SW.EC.1..,E16)8:,31 

30  TT=X?min+{NUP*0ELX2> 

IF(L.EQ.l) TT=TT-0ELX3 
PA0=S*NHCF*CSLX1*EXP(ZXP(TT> > 

GO  TO  1 

31  RAD=NH0R*0ELR*l.CS-:5 

1  IF(MPASS.C-T.N$TCF)GO  TO  15 
K=MOCE (NPASS) 

M=LMAF 

IF(LMAP.£Q.3)M=1 

CONSTRUCT  CONTOUR  CURVES  FROM  PCINT  VALUE  RATA. 

A=VUL(NoasS) 

KK=1 

2  CALL  CONTOL::,(m,K,i<K,NuOp,NUP,A,S,CELX?,CELX1,  X1MIN,  X3MIN,  X3SW ,  DELF 
A,NN,rATft ,N,F) 

SET  U?  THE  PLOT. 

CALL  SETUP (”, K, RAD, SCALER, SC ALEZ, WORK A) 

HI=ALT*SCALEZ 

DRAW  THE  PLOT 

LIMIT=2 

IF(M.EG.2) LIMIT  =  16 
LCOUNT  =  j 
X  =  3.  8 
Y=  .  3 

CALL  S  Y  M  3  0  L  ( j  •  0  » H I  , g . 1, 11 ,0 . C ,-l) 

CALL  S V M 3 O L ( X  ,Y  , j . 1, 11 , C . : , - 1) 

CALL  SYM3CL (995. ,999. ,  :.1,6H  °URST,C.A,6) 


125 


srsrwaw 


'fltwnwr 


1 


GNIi/PH/72-8 


4 

10 


12 

11 


20 


5 

25 


27 

28 

f> 

q 


7 

15 

lfi 


LC0UNT=LCCUNT+2 

IF (LCCUNT .GT.LIMIT)GC  TO  5  ' 

J=LC0UMT/-2 

NP0=M ( J) 

IF(NPO.EG.l) GC  TC  2 
00  U  1=2, NFQ 

C=QATA(I,LCC'JNT-1)*SCALEP 
0=0ATA  ( I  ,LCC'JNT)  'SCALE7 
CALL  SYMBOL (C,0 , 0 .1, J, L. b ,-i> 
Y=Y+ • 2 

8=0 AT  A(i,LCOUNT-l) 

IF  (M.  A'c.  1) GO  TO  12 
90=ALOG13 ( 2) 

NB=3P 

AP=1C 

3=B/(98**N«) 

GO  TO  11 
NN=NN-1 


3N=NN 

BA=1£ 

CALL  SYM-30L(X,Y,3.1,J,C.t,-i) 

IF (M. EO. „ ) GO  TO  2 C 

CALL  NUMt’lIR(X4.2,99Q. ' ,G.i,9A,C .0,-1) 
CALL  NUMBER (999. ,Y+U. . 5, 3 . 1, BN, o . ,-l ) 
GO  TO  3 

CALL  NUMBER (X4. 2, 999. 3,0.1,8,0.1,41) 
CALLSYM80L(X+.7,Y  ,  .1, 4,  j  .  ,-l) 

CALL  NUMBER (999. ,999. , u. 1, A3, T . ,-l) 
CALL  NUMBER (999., Y4C. -5, 6.1,03, C.,-1) 
GO  TO  3 

IF(M.fQ.2)SC  TO  28 
NCOUNT  =  3 
NC0UNT=NC0UNT+1 
IF(NC0UNT.GT.NAC  )GO  TO  27 
IF (A  PRAY (NCOUNT) .EC. J ) GO  TO  25 
IF(P0P(1, NCOUNT)  .5T»2-DGr  TO  25 
I=NC0UNT  +  5  4 


C=P0S(t, NCOUNT) ^SCALE3 

0  =  P0S(2,NCCUNT) *SCAL£7  .  .  , 

CALL  SYMBOL ( C, 0 , 3 . 1 , I , K . l , -1)  *■  i 

GO  TO  29  1 

CALL  SYM30L(2.9,3.1,'3.1,19H1,2,3,ETC.  AIRCRAFT, C .u , 19) 

CALL  PLOT (12. .,-2. C, -3) 

IF(LMAP. EC, 7) 6, 7 
IF(M. £0. 2) 7,9 


M  =  2 

KK=2 

GC  TO  2 

NFASS=NPA3S41 

GC  TO  1 

NCOUNT  =  r> 

NC0UNT=NC0UNT41 

IF (NCOUNT. GT , MAC  )G0  TO  5P 

IF{AF°AY( NCOUNT) ,tC.4) AOPAY(NOCUNT)=1 

IF (A»PAY (NCOUNT) .EO. 2. OR. ARRAY (NCOUNT ) . EQ 
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GO  TO  16 
REWINn  n 
RETURN 
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Subroutine  CONTOUR 

This  subroutine  calculates  points  on  isofluence  or 
isodose  lines.  The  subroutine  is  constructed  so  that  the 
entire  meshed  area  is  scanned;  however,  if  any  points  appear 
that  will  be  oux  of  the  range  of  the  plot,  this  point  will 
not  be  stored. 

Subroutine  CONTOUR  Glossary: 

A 

DATA 

F 

NPO 


The  fluence  or  dose  value  for  which  the  line 
is  being  calculated 

The  results  of  the  mesh  scan.  This  array 
contains  the  r,z  coordinates  of  the  points 

The  total  fluence  or  dose 

The  number  of  points  defining  an  isofluence 
or  isodose  line 
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SUBROUTINE  CCNTCUF  <M,!<,t<KtNHCR,NUP,  VLL,S,0ELX3,CELX1,X1MIN,X3M 
AX?SW,rELR,»CN,raTfl,K,F) 

rCMMOfi  a,l.IMIT,B,II,J,KO0,Kl,Jl,I,C,X3,Z,R,0EL,Xi 
DIMENSION  DATA  ( 2  C  B  ,  16)  ,N(fi>  ,F  (6(3,120 

THIS  Sl;P"CUTT\E  SCANS  THE  MESHEC  AREA 
AMD  DETERMINES  ,FOR  ANY  GIVEN  VALLE  OF 
FLUENCE  CF  rcSE  ,WHERE  THE  VALUES  EXIST  IN 
TERMS  CF  RADIUS  AND  HEIGHT 

if(kk;ne.dgc  to  i 

^EAB  ( ) (  (F  (I,J) ,I=l,NHrF), J=i,NUF) 

1  GO  TO  (2,’)M 

2  A=VUL 
J=G 

LIMTT=2 
GO  TO  3 

3  LTMTT=16 
3=F(i,i> 

00  A  1=1, NLP 

4  B=AMAXl(9fF(l*I) ) 

NN=Q 

5  QrR/1* 

NN=NN+1 

IF (B*  LT« 1, ) 6  >  5 

6  ‘  NN=NN-2 

A=10 « **NN 
J=0 

7  A=A/1C 

8  NPO=i 
J=J  +  2 

IF(J.GT.LIMIT)GO  TC  40 
Kl=NUF-l 

DATA  (NF0,J-1)*A 

NF0=NF0+1 

DC  9  11=1, NHCR 

DO  1C  J1=2,K1 

IF(F(II,J1) .IT. A) GO  TO  1C 

TF(F(II,Jl-i> .GT. A) GC  TO  11 

0EL  =  pFLX3*(A-F(II,Jl-l))/(F(IIf  Jl)-Fni»Jl-l)  > 

X?=X3MIN+(CELX3*  < Jl-2)  )  +OEL 
GO  TO  12 

11  IFCFdl,  Jl«-1>  .GT.  A)GO  Tn  in 

0El=CrLX3*  '  A-F  (IT .  J 1 ) } /(F (II,J1+1)-F (II, Jl) ) 

X?=Y?MINMCElV3*  (Jl-1)  )  +OEL 

12  ?=S*EX°(X7) 

IF(X3.GE.X3FW)GC  TC  13 
X1=X1MIN* (CELX1* ( I I - 1 ) ) 

R=S*Xi*EX? (EXP(Xl) ) 

GC  T0  14 

1?  R=(XlMIN+( ( 1 1-1 ) *  CELR ) )  *1  •  GE-OF 
14  IF (R , GT , 2. 0 . >  GC  TC  10 
DATA (NFC, J-1)=R 
DATA (NPB,  J)  =  ”> 

NFC=NF0+1 
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IF(NPO.EG.2CO>GC  70  21 
ln  CONTINUE 
«  CONTINUE 
X1=NH0P-1 
00  15  J1=1,NCP 
X3=X3MIN*(C£LX3* ( Jl-i)  ) 

C-0ELX1 

IF(X3.GS.X2SW>G=CELR 

7=S*EXP(X3) 

or  If  I=2,K1 

IF(F(I,J1> .LT.a)CC  TC  16 
IF(F(T-1,J1)  .C-T.fl)GC  TC  17 
0EL=D*(A-F(I-1,J1))/(F(I,J1)-F(I-1,J1)) 
X1=X1MIN+(C*  (1-2)  WCEL 

17  IF(F(T  +  1,J1)  .GT.  A)C-C  TC  16 
OEL=0*(A-F(I,J1))/(F(I-H,J1)-F(I,J1)  ) 
Xi=XiMIN* <G*  (1-1) ) +CEL 

2G  TF(X3.GE«X2SW)GC  TC  18 
R=S*X1*FXF(EXP(XT)) 

GO  TO  ig 

18  °=Xi*l • QE- T5 

19  IF(R.GT.2? T.)PC  TC  16 
OATA (NF0,J-1)=P 

DATA (NP0,J)=7 
NFO=Nc  J  +  l 

IF(NPC.EC.2Cj)GC  TO  21 
16  CONTINUE 
IB  CONTINUE 
21  J1=J/? 

N(J1)=NP0-1 
GO  TC  7 
*»n  RFTUPN 
END 
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Subroutine  SETUP 

This  subroutine  completely  sets  up  the  plot.  The  axes 
are  drawn  and  the  correct  axis  labels  are  selected  and 
drawn.  In  addition,  this  subroutine  selects  the  plot  title 
based  on  what  type  of  plot  is  to  be  drawn.  The  scale 
factors  required  by  MAP  for  plotting  are  also  calculated. 
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SUBROUTINE  SETUP (H, MODE, RAD, SCALER, SCALE!, KORKA) 

DIMENSION  HORffl Cl "’241 

COMMON  <C, MM,  L, N,  Y,X,0,CA,  I, 0,E, DO 

C0MM0N/TITLES/IA'11,3),TC5,22) 

THIS  SUBROUTINE  DRAfcS  AND  LABELS  THE  AXIS 
ANO  TITLE?  THE  FIGURE. 


FIRST  LOCATE  THE  ORIGIN. 

CALL  PLOT <2. :,-7. ^,-3) 

CALL  PLOTC j. 0,1. 5,-3) 

CALL  PLOT C  6. j ,G. 1,2) 

CALL  FLQT(6.j,3.45,2) 

CALL  PLOTC*. 1,3.45,2) 

CALL  PLOTC  0.3,fJ.a,2) 

'  SET  UD  THE  FIGUFE  TITLES. 
K=MCOE 

CO  TO  C1,2)M 

GO  TOf?, 4, 5, 5, 6, 5, 5*6) MODE 

MM=17 

GO  TO  9 

MM=18 

GO  TO  9 

MM=19 

GO  TO  9 

MMt2C 

GO  TO  9 

K=MODE+8 

GO  TO  (7, 7, 7, 7, 3, 7,7,3)MOCE 

MM=21 

GO  TO  9 

MM=22 

SET  UP  THE  AXIS  LABELS. 


0=6HHE IGHT 
E=6HRA0IUS 
00=5H  (KM) 

OETER( INE  THE  RADIUS  AND  HEIGHT  MAXIMUM 
COORDINATES •  DETERMINE  THE  RADIUS  AND 
HEIGHT  scale  factors. 

IF(RAD.GT.5:.)GC  TO  1 ' 

N=5 
L  =  1 

GO  TO  12 

IF{°AD.GT.iu  ' « )  GO  TO  11 
N  =  1C 
L  =  2 

GO  TO  12 
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11  N=2C 
L=3 

12  SCALE7=5./i:i!. 

SCAL£R=5 •/  (K»13.) 

DRAW  THE  AYTS. 

CALL  PLL’T(  J.9,~.42,-3) 

CALL  PLC'( 3. J,r.n,2) 

CALL  PLOT( j.o, j. 3,3) 

CALL  PLCT(5. :9u,1f?) 

CALL  PLCT< 3.J, 3.3,3) 

DRAW  THE  TIC  MAPKS. 

Q=8. /I J  * 

Y=0. 

DO  13  1=1, il 
CALL  PL0T(„.,Y,3) 

CALL  PLOT ( -3 • j7 , Y  ,2 ) 

13  Y=Y+0 
GA=5./1G. 

X=G. 

00  14  1=1, 11 
CALL  PLOT (X,j.,3) 

CALL  PL0T(X,-C. 37,2) 

14  X=X+Q A 

NUMBER  TH£  AXIS. 

Y=-« . 1 

00  15  1=1,11 

CALL  SYMBOL (-. 4, Y, L.l, IA  (1,2) ,C . ,3) 

15  Y=Y+0 
X=-.2 

00  If,  1  =  1,11 

CALL  SYMR0L(X,-.2,C.1,IA  (I,L)»r:.,3) 
If  X=X+0 A 

LABEL  THE  AXIS. 


CALL  SYH30L(i.75,“.4,' .13,E,C.  ,6) 

CALL  SYMBOL(939.,99Q.,C.i3,DQ,*.0,5) 

CALL  SYMBOL (-.5, 3.25, 3 .13,0, 9.  . C, 6) 

CALL  SYMBOL (399., 999. ,«.l?,00,«u.w,5) 

LABEL  THE  FIGURE. 

CALL  SY^f30L(--.9,-C.f5,u.l3,T(l,'<),L..,46) 
CALL  SYu*0L(-'.9,-C.95,*.13,T(l,MM)  ,0,C,45) 
RETURN 
ENO 
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BLOCK  DATA 

This  subroutine  stores  data  into  labeled  common. 

BLOCK  DATA  Glossary: 

T  An  array  that  contains  the  titles  to  all  the 

plots  that  can  be  drawn  by  subroutine  MAP 

IA  An  array  containing  three  sets  of  coordinates 

for  the  plot  axis 

SPECFN  The  fission  weapon  neutron  output  spectrum 

SPECFG  The  fission  weapon  gamma  output  spectrum 

SPECTNN  The  thermonuclear  weapon  neutron  output 

spectrum 


SPECTNG 


The  thermonuclear  weapon  gamma  output  spectrum 
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9L0CK  DATA 

COMMOM/S°ECTRA/SFECTNN  (22)  »SPECFN(22)  ,SP£CFG(18)  ,SPECTNG(1B> 

COMMON /TIT  LES/IA (11 »  3)  ,T<5,22) 

DATA  <  (T<J,J)  ,1  =  1,5),  J=i»10)/lv.H  FIG.  ,10P  NEUTRCN  , 

UICHFLUEN'CE  VU,1 JHLNERASILIT, 

ASHY  ,  1 JH  FIG.  ,1GH  GAMMA  F, 1GHLUENCE  VUL,  nFNERAnILITY, 

B6H  ,  1  OH  FIG.  ,iCHNEUTRON  TI'9iqhSSUc  COSE  »19FVL'LNERA°IL, 

CGHITY  ,10H  FIG.  ,10H  GAMMA  TIS, 10HSUE  DCSE  V, 10HULNERAEILI, 

06HTY  ,i3H  FIG.  ,1CH  NEUTRON, 1GH  +  GAMMA  T,1CHSSUE  DCSE, 

E6H  •  ,1 jH  FIG.  ,1CHNCUTR0N  SI,1CHLIC0N  COSE, IGF  VULNERA°I, 

F5HLITY  ,10H  FIG.  , iCH  GAMMA  SIL,ltHICCN  COSE  ,  19FVULNE0AEIL,i 
G6HITY  ,  1  OH  FIG.  ,i»*h  NEUTRON, 10H  +  G®  MMA  S,10FILICON  COS,i 

H6H5  ,13H  FIG.  ,1CF  NEUTRCN  I ,1CHSCFLUENCE  ,13FLINES  (N/C,> 

I6HM2)  ,  1  JH  FIG.  » 1GH  GAMMA  IS , 1GHOFLUENCE  L,10HINES  (G/CM,! 
J6H2)  / 

DATA  <(T(I,J> ,1=1,5), J=11,22)/1CH  FIG.  N,1GHEUTR0N  TIS, 

X1CHSUE  IS0003,10HE  LIMES  (R, 

K6HA0S)  ,  1  JH  FIG.  ,1CFGAWA  TISS,luHUE  ISCDOSE,i3F  LINES  (RA, 

L6HDS)  ,  1  OH  FIG.  ,10F  NEUTRON  ,1C-H+  GAMMA  TI,li3HSSLE  ISCDC, 
M6hSE  ,ijH  FIG.  N,lrHEUTROM  SIL,UMICCN  ISODO,1GFSE  LINES  (,, 

Y6HRA0S)  ,  1  OH  FIG.  , lrFGAMMA  SILI,ltHCON  ISODCS,10FE  LtNES  (R,, 

W6HADS)  ,  1  OH  FIG.  ,1!F  NEUTRCN  4,ltH  GAMMA  SIL,1JFICCN  ISCOC,i 

N6HSE  ,13H  L, 1CHINE  <N/CM2,1GH)  EXPONENT , 13FIAL  AIR  , 

06H  , 1 GH  L,1lFINE  (G/CM2,1GH)  EXPONENT, 1C  HI AL  AIR  , 

P6H  ,  1 QH  LjiOHNE  (RA0S),1(jH  EXPONFNTI ,  1C  HAL  AIR  , 

OSH  , 1  OH  VULNFRA9,1CHILITY  LINE , 1GH  (RADS)  EX ,10 FFCNENTI AL  , 

R6HAIR  , 1GH  ,UF  EXFON, 1DHENTIAL  AIR,10F  , 

S6H  ,1GH  L,lr>IN£S  (RADS, ICH)  EXPONENT ,  10  HI  A  L  AIR  , 

T6H  / 

DATA  IA/Ori  G,3H  5,3«  16, 3H  15, 3H  2J,3H  25, 3H  3w,3H  35, 2H  AO, 

A3h  45, 3H  5j,3H  0,3H  lu,3«  23, 3H  30, 3H  40, 3H  5C,3H  6J  ,"*H  7j,3F  80! 

93H  9C,3HlvQ,3H  9,3H  ?C,3H  40, 3H  60, 3H  8 j , 3H1 or , 3H12 J ,3H14 C , 3H16 0 
C3H18C , 3H2C  0/ 

DATA'  SPECPN/3.9  2EU9,2.232E  +  2J  ,  8 . 7E  + 2  C  ,3. 48E  + 21  ,  8 . 7^.5  E  +  C  1 , 8 . 7 '‘5E+ 
A21, 3*1. >+951^22, 2*4, 23E  +  22, 2*4, 2325F  +  22,3.875E  +  21,3*C.'’/ 

DATA  SP5CTNN/6.0Jlt+22,2.176E+c?,1.1985E+22,1.2495E+22, 
A2*1.4675E+22, 3*1.41676+22, 2*3.625E+22, 2*7. 9475Eh22,?.1C?5E+23» 
98*0.0/ 

DATA  SPECcG/i.',54  8t  +  19,3.,:]i9-  +  lS,1.3539E  +  2w  ,  4. 7555E  +  2L-  ,4. 75  5<:E  +  2 
A,1.C71SE+21,1. j71C£+21,?. 3562E+21, 3. 22E+2i,4. 08385+21,3.7741 E +21,, 
34. 512l+21,5.2A99E+21, 2*2. 2981F+21, 2.7 '62E+21, 2*0,0/ 

» 

OATA  SPcCTNG/G.  32^93  +  18, 2, 9F0Pr+19, 5. 2693E+19, 2*2. 377s"  +  2rl, 
A2*5.3595E+2C»1.17alE+2i,1.615r+21,2.J419E+21,t.fi37E+21,2.255°E+21' 
92.6249F+21,2*i.i49F+21,1.3531E+21,2*0.C/  S 

end  ; 
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Appendix  C 

Sea  Level  AiT  Cross  Sections 

The  following  pages  contain  the  40  group  sea  level  air 
cross  sections  used  in  the  code.  This  page  gives  the 
instructions  on  how  to  read  these  pages. 

The  data  is  listed  in  40  columns,  one  for  each  group. 

The  first  column  is  neutron  group  one,  the  second  column  is 
neutron  group  two,  and  so  on.  Since  the  neutron  and  gamma 
cross  sections  are  coupled,  the  data  can  be  regarded  as  one 
40  group  cross  section  set.  Therefore,  column  23,  which  is 
gamma  group  one,  can  be  considered  as  group  23.  Each  column 
contains  43  values  for  group  cross  sections.  The  first 
number  is  the  group  transport  cross  section,  the  second 
number  is  the  group  removal  cross  section,  the  third  number 
is  the  group  total  cross  section,  and  the  fourth  number  is 
the  in  group  scatter  cross  section.  The  remaining  numbers 
in  the  column  are  values  for  the  group  to  group  scatter  cross 
sections.  So,  the  fifth  number  is  scatter  to  that  group  G 
from  group  G-l,  the  sixth  number  is  scatter  to  that  group  G 
from  group  G-2,  and  so  on.  All  the  cross  sections  are 
macroscopic  with  units  of  cm  *. 
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